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DIFFERENCE  METHOD  IN  FLUID  DYNAMICS 


I.  NUMERICAL  SOLUTION  OF  TWO-DIMENSIONAL 
VORTICITY  EQUATION 

Guo  Benyu 

Shanghai  University  of  Science  and  Technology 


ABSTRACT 

With  the  two-dimensional  vorticity  equation  as  an  example, 
this  paper  systematically  discussed  the  theoretical  problem  of 
difference  equations  in  fluid  dynamics.  Many  schemes  are 
constructed  In  this  paper  based  on  the  conservation  laws  and 
transport  properties.  Generalized  stability  is  introduced. 

The  error  of  periodic  solutions  is  rigorously  estimated.  The 
effect  of  error  in  boundary  values  on  the  stability  of  compu¬ 
tation  and  the  boundary  shape,  boundary  conditions  and  its 
treatment  in  the  initial/boundary  problems  are  analyzed. 

The  computational  problem  of  non-viscous  flow  and  large 
gradient  solution  is  treated  according  to  the  conservation 
laws.  Lastly,  the  convergence  of  dynamic  relaxation  and  the 
existence  of  solutions  of  steady  flow  computational  schemes 
are  studied.  All  the  results  in  this  paper  are  proved  rigor¬ 
ously,  and  the  method  used  is  also  applicable  to  the  general 
fluid  dynamic  problems.  Examples  are  used  to  illustrate  the 
applications  of  these  results  In  fluid  dynamic  computations 
and  numerical  weather  forecast. 
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The  theoretical  study  of  difference  method  basically 
consists  of  five  aspects:  (1)  the  principles  of  scheme; 

(2)  stability  of  solution  against  perturbation  on  the  RHS  and 
the  initial  values;  (3)  non-viscous  flow  calculation; 

(10  effect  of  boundary  shape,  type  of  boundary  condition,  and 
error  in  boundary  value  on  the  solution  of  initial/boundary 
value  problems;  (5)  steady  flow  calculation  [1-5].  A  large 
amount  of  work  has  been  done  in  foreign  countries  but  there 
are  only  a  few  very  special  theoretical  results.  Based  on 
the  work  of  Ref.  [6-8],  we  discuss  in  this  paper  all  aspects 
of  the  five  problems  quoted  above  with  the  two-dimensional 
vorticity  equation  as  an  example,  and  obtain  a  more  systematic 
result.  The  method  is  applicable  to  the  general  fluid  dynam¬ 
ical  problems. 


I.  Symbols  and  Lemmas 

R  represents  an  open  region  in  (x-^,  x2)  plane  with 
boundary  r.  The  grid  step  length  along  Xj  direction  is  h. 

The  coordinates  of  the  grid  point  Q  are  *,tp)  —  *(P)*>  where 
kj(Q)  is  an  integer.  The  coordinate  of  Jfc.CP)  is 
*,(P*'0  -  UP)  ±  *....  represents  the  set  of  internal  grid 

points  with  boundary  r^.  (or  Tjm)  is  the  set  of  boundary 

points  that  makes  Q-*i  (or  Q*‘i  )  belong  to  R^.  RjM  (or  Rjm) 
represents  the  set  of  all  internal  points  at  distance  h  from 
TjM  (or  Tjm) .  T,  —  r,M  4-  r,„ ,  R,  -  +  J?,„  R J  -  /?,  +  R,.  t  is  the 

step  length  along  t,  l  —  r*"\  «'(P.  O  represents  the  value  of  the 
net  function  w  at  time  kx  and  position  Q.  Sometimes  it  is 
simply  denoted  by  wiK)  or  «>.  w  —  *?($)  —  «k(^  +  1). 

represents  the  forward,  backward  and  central  difference  coeffi¬ 
cients  along  *,  .  w.{Q)  is  the  outer  normal  difference  coeffi¬ 
cient  on  .  When  P€R,«  ,  *•.<  p )  ~  »,/(.(> )  .  When  QtR..  , 

».( P)~  —  *'«,(P).  represents  the  tangential  difference  coefficient 
on  R*  *t  .  When  p  t  Rim  h  Tn,  i  a'i(O)  ■■  *t'i,(p)»  •  When  q  €  R +  rtm  , 
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“'•(£>)"  —  u>,(Q)  .  When 


«'i(CO— —  When  Q  (  R>m  +  , 

p  €  Rlm  4-  r,„  ,  «-,(p)  —  w'i^P) 

If  PCPJ  ,  we  shall  denote  by  the  point  on  f*  at 
distance  h  from  Q,  and  also  denote  £  1«<0'MP)  +  «'(£>)«•({?’)] 

simply  as  X' (“•*'  +  “'«' )  etc. 

We  further  define 

l  1  1 

A*.  *«■■“-(  +  -  (viv,,)..,  A'tt/  -  5]  A' .«/ 

*  "  •  m  I 

Denote  „*•,  —  p**  >,  is  an  ADE  type  operator,  i.e. 

H*( «*).—  WX«)  is  taken  according  to  a  certain  order,  where 


—  —  i  (v  +  v"**V.i,  —  -7 +  V'Ow.t, 

H)C  iv)  —  Y  (»  4-  v+*< —  y  Cv  +  v",0«»i«l 

«;(«-)  —  —  (v  +  +  ”•(*’  + 

2  2 

HRu>)  —  —  (y  +  ~  (y  +  yt,,)«'i*,:; 

2  2 


We  define  the  following  internal 

C«y.  y)*“  A'^CpMp), 
»•«* 


products  and  norms: 
M*  “*  (w,  w ) 


;  S<llv^y'.,|!‘ +  llVrv*«yi<l!,>  ,  if  v  »  1  ,  then  simply 

•  (•! 

denote  by  fln-BJ 


ll-  Ji-T  £iOkiK  +  II 

«  (-1 


"\\U  - 

»*r, 
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Denote  also 

**  M'l  “  —  («  U  +  uv'')w,  C„*C“|  v»  w)  “  yi  ®*/C U;  r •  **') 

'  2  V  * 

v 

>  2 


D.W)  -  S  (-7  «'**  -  7  «^) 

*i«f*i*  2  ... 

—  2  (—  k'1  rV*''  —  —  w'1) 

«,.***.' 2  2  ' 

—  |MI!  +■  r  +  y)  (|!«'|l!,)i  ~  0r|l“'»lll  +  frlMlo 


+  MU  +  +  mIM* 


where 


_  «  _  ff  +  *£L, 
2  2 


l'»U)fc.  -  + 1 S 

/,. «.  o’-  |.*»r  +  r  2  cuxor  +  uxor  MMW 

/  •*  0 


It  can  be  proved  that 


2(«r,  «/,)  -  Ctl«'I!*>»  “  tlM* 

(1.1) 

2(»r>  u/|)  (ll«'ll,)i  +  f IIm'iH* 

(1.2)  /131 

Kw,  »„><-  imij,  -  n^.ir  —  i».r 

(1.3) 

(1.4) 

?f  •  A*f  <  f  (»«’.,♦  u.,)  +  Cv»i;,  r.,)  —  2Bk.Cv,  O  —  0 

(w,  A’«0  +  flu'll*.  —  B*;(r,  «u,  w.)  —  0 

(1.5) 

.  » 

2{*i,  A’*/)  +  (Bu  ll?.),  -  tDw.B!.  -  —  2j  (»’•»  +  “'•/> 

<1-1 

—  2B«j(u,  «r„  IT.)  —  0 

(1.6) 

+  (BwlH.)*  -  4"  22  <*•  “'i/  + 

2  |*t 

—  iBtfiv,  »,  «r„)  •-  0 
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(1.7) 


i<*„  «'(«/))  +  AIU.II?,  -  (»i,,  «-,>)  -  -  0 


(1.8) 


Lemma  1.  If  ?  >  ^  >  0,  >  then  J!«'*!l*  < 

Lemma  2.  If  Hr,  —  fir,’  -  0  ,  then  ||«,<,||>  AIIii'IIMMUUIIl 

Lemma  3.  If  Hr*  —  o  ,  then  IMIJ <  HAwl!1 
Lemma  4.  Assume:  (1)  «•(*),  KO  is  a  non-negative  net 

function.  p(T)  is  non-negative,  t  is  suitable  small, 

Mu  v  »:<  «-,«•  •  *My  +  tt'(/OA--)+/*l«’(A))‘'(0  f 

,  /•(*')  <0;  (3)  «*-(0  )<•-#>(  7").  when  P(T)  +  r  S  /!-(»)*.  QO 

1  *■  J  , -0 


pr**"*  r.‘  fh'  .  ,  then  when  *r  *S  T  *£=  T,  ,  «-(0  ^ 


II.  Basic  Principles  of  Scheme  Construction 


Let  {(*!»*•»  0  and  M*n  *«,  t  represent  vorticity  and 
flow  function.  Then  the  2  dimensional  vorticity  equation  is 


S-£2:+ £  £-&£(•£)-*• 


VV  -  f  -  /„ 


X  |0,  T] 


(2.1) 


Where  <)  is  known,  *«,  O  is  coefficient  of  vis¬ 
cosity,  <  We  ais0  denote 

K*  “*  inf  KP),  #*  —  mp  ^'(P)>'_I(P,) 

p,** 

and  assume  |~|  and  |— l|  are  bounded. 


The  principle  of  scheme  construction  is  to  express  physical 
laws  with  suitable  separation  form,  for  example,  if  »>*=/,  *^o  , 

then  we  get  from  (2.1)  the  conservation  laws  for  vorticity  and 
its  square: 

*  *f.»*  "  Jj,  *•’  *i. 


/132 


'  +  ir/")'r  '-1*2 

As  for  the  difference  scheme,  since  the  step  length  is  non¬ 
zero,  we  oannot  simultaneously  simulate  both  conservation 
laws.  The  K,  conservation  law  is  simulated  by  the  Lax  scheme. 

r  r  s  i  ^ 

V.'eL  J  simulate  the  conservation  property.  Let  us  define 

jl  <  r,  Uf)  •*  »  »•(,,  Ja,‘ l(*'>  **  •'I,*4'*, 

j:  ,c». « )  -  C4'*4'*,).,.  /:,.»(«'» “')  •=  («"4'i,>», 

*4')  -  -(•'i,44').,.  «'>  “ 

«')  — J*. 

j,\  it*.  <*)  “  -(*'«4'i,)»,»  •/!,•»<> »  *')  “  -(*'*4'<>. 

«*)  -  c*'!,44').,,  /;,•»(«',  *0  -  Ui.nOi, 

7.r.<*.  «0  -  y  «0  +  y  “') 

*  ^ 

y/'*,  •-)  -  2  *0»  /7C».  <0  “  £  *0 

i-i  <•> 

jt c«<,  k') — 4"  *«') +  4" 

2  • 

»  » 

re*. «-)  -  2  »).  J"(*»  «0  -  2  «*> 

*•» 

•  * 

y.(y,  -  >  -  2  *<),  yu.  v)  -  2  y*/-*'* 

»i  !•> 

where  O(^o,  2  ®<  *"  1  *  ***  can  Proveci  that 

*•  i 

*•*»  y., ■>(*-.  •<  ))  +  C*.  j.rt(u,  mO)  -•  V ,  1) 

(*.  y.t*. « ))  +  (•>,  «/))  «■  — b*«(«,  u'Pt,  i)—  b *•(•>,  i) 

If  «.  —  a,  ,  then 

C*,  f(  p,  *■))  t  (»,  J(u,  |*»))  —•  a,B*»(«li';,  P,  1)  +  Oj6vi*( I'M*,,  «,  1 ) 

~  •»**•<•».  1)  —  a.B«*(r,  vmi,  l) 


(2.4) 


i 

» 

i 

i 

« 

i 

i 


(2.2) 

(2.3) 
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If  one  the  the  2  following  conditions  are  satisfied  on  the 
boundary  /'»  : 

A.  The  various  quantities  are  periodic  with  periodicity  L„  Rt 

along  the  direction  .  R^  is  the  rectangle  (0,  L,  —  A;  o,  L,  —  A); 

B.  is  a  limited  region  and  the  flow  along  1\  is  zero. 

For  example,  u  —  v  —  w  —  0  ,  then  the  RHS  of  equation  (2.M  is 

zero.  If  we  now  use  i<  <P  and  A*.  <j>)  to  approximate  £»  &  and 

the  Jacobi  operator  in  equation  (2.1),  then  when  a*  **  a>  , 
the  conservation  law 

(.ij.Ai/.  <p))  —  >2,1)  —  <pn»  l)  ^ 

is  well  described. 

If  condition  A  or  B  is  satisfied  on  the  boundary,  and 
,  then 

(^,  7,/ij,  <p))  ~  0  (2.6) 

The  weighted  average  conservative  scheme  for  computing 
equation  (2.1)  is 

O'n.t  J{ n  +  8th,,  ?)  -'<&’(  ij  +  0T1{ —  Xrh~'H,{ri')  “  /i  1 
L,(i.  <p)  —  by  -  h  —  1,  1(2.7)  /133 

where  e  <o,  <  a,  1,  o  $(0  <  i,  o<  Z(0  <  1,  —  -1  <  eCO  <  o  , 

but  when  'M  >  is  calculated,  we  always  take  SCO)-0  .  In 
the  explanation  below,  we  might  as  well  assume  that  condition 
A  is  satisfied. 

1.  The  function  of  5  is  to  maintain  the  conservative 
property.  If  fact,  if  *  ■«  /.  —  0  —  *  —  •  —  0,  a,  —  «,  ,  we  still 

have  IliHJ  —  Tliij.ll1  .  Hence  each  step  of  the  computation  will 
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see  a  small  increase  of  the  virtual  energy  which  will  cause 
energy  explosion  and  computational  overflow  when  it  exceeds 
some  critical  value  through  long  term  accumulation.  But  if 
8  “  1/2  ,  then  M*  "* 0  and  the  conservation  is  strict. 

2.  The  function  of  0  is  to  filter.  It  equivalently 
filters  the  n  in  to  be  V*  “  *}  —  .  if  the  harmonic 

function  of  the  computer  error  >1(0  “  Ac‘rk  ,  then 

fl’CO  -  (l  +  40sin ‘^JiKO 

Because  —  -  <  6  <  o  }  hence  l|»J*l)  <  M  ,  so  that  its  growth 
is  suppressed.  Robert  has  used  explicit  filtering  and  0  "0.01 

See  [14]. 


3.  If  X  *  0  ,  then  equation  (2.7)  is  the  successive 

overrelaxation  scheme  which  is  used  to  improve  computational 
stability.  Usually  the  WftO  are  used  in  turn  to  obtain  fl<> 
and  then  let 

-  7  X)  * 

which  is  the  Larkin  method^^’ 

4.  In  (2.7)  we  can  also  use  the  linear  combination  of 

f^  and  f^  to  replace  fi*  .  All  the  results  in  this  paper  will 
also  hold. 

The  weighted  average  conservation  method  may  be  used  In 
conjunction  with  the  splitting  method,  namely  with  n*  as  the 
supplement  value. 

a;t(*  +  ?*)  +  j  U 

*  *  ...  , 

«  -  ^  + 1, »)  +  y +  f)  +  j /,  (2>8j 
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If  /,  b  t.  i=  o, «.  »  a,  -  -  ,  then  from  equation  (2.6),  11*711*  ■*  U»;ll' 

Since  equation  (2.8)  not  only  strictly  satisfies  the  conser¬ 
vation  laws,  but  also  may  be  computed  explicitly,  it  is  to 
be  recommended. 


Another  principle  for  scheme  construction  is  the  trans¬ 
port  property.  Let  Q  be  an  internal  net  point.  Form  tne 
triangle  with  £>,  as  the  vertices,  with 

Q,  Q*'‘.  Q'‘‘  as  vertices,  and  similarly  for  £>„  Z),.  .  If 

<p«,  <0,  <p,(<0.  ,  then  the  wind  blows  toward  Q  in  Di  .  Denote 
it  by  «(D„  v>)  —  1  else  <p)  >-  0  .  Similarly  the  other  £(£>,,</>) 

may  be  computed.  We  shall  construct  the  difference  operator 
F0j»  <pj  f  or  the  following  three  cases: 

1.  e(L>,,  y;  —  s,,'  ,  indicating  that  the  wind  blows  uni¬ 
directional  ly  toward  Q.  For  example,  if  /a"*2-  ,  then 

F(y>  <p)  — 

2.  There  are  several  £(D,.'c)“'  ,  indicating  that  Q 

lies  between  several  opposite  currents.  For  example, 
e(D,,  <j>)  —  »,.i  +  ,  then  let  F(v,  <p)  —  rF  fa,  tp)  +  O  —  r)F,(ij,  <p) 

in  which  o<r<l,  F.C^,  <p)  -  <p.,n.,  -  V’,1)'.-  . 

3.  All  e(D(,  9>)  —  o  ,  indicating  that  Q  is  the  center  /13^ 

or  source  of  vorticity,  and  therefore  F(>j,  <p)  —  o  .  From  this 

[7  22] 

may  be  obtained  the  modified  counter-wind  scheme  ’ 

tp)  —  t),  +  Ott),,  —  F(.jj,  i p)  —  ^*,C»7  +  any,)  —  «=  /,  l  (2.9) 

Li(v>  <p)  “  f)  » 

III  Error  Estimation  in  Problems  with  Periodic  Solution  and  the 
Generalized  Stability 

Let  the  difference  scheme  be  —  /(O  where  /(*) 

is  the  condition  for  solution  determination.  The  error  of  ] , 
will  Induce  an  error,  >7,  in  v  .  Ordinarily,  stability  means 
that  there  exists  an  absolute  constant  M  such  that  ll*/ll  <  Mll?ll 
However  most  non-linear  schemes  do  not  have  such  properties. 
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f  7  21  221 

The  author  ’  *  suggested  a  generalized  stability, 

namely  that  there  exists  positive  numbers  M,  N  and  constant 
«'  independent  of  h,  r  such  that  when  Nh‘\  kr  <  T  <  T.(!!?il)  , 

||q||  A/ ll?il .  The  highest  lower  bound  is  called  the  index  s  ol‘  the 

T  2 1 1 

generalized  stability.  Some  explanations  are  offered  belowL 

1.  For  non-linear  mechanical  systems,  generally  only 
when  H/ll  satisfies  certain  conditions  will  there  exist  a 
unique  solution  which  is  only  stable  against  perturbations 
within  certain  range,  e.g.  the  3  dimensional  Navier-Stokes 
Problem.  But  if  L»o  —  0,  <  <  o  ,  then  when  ll/ll  <  N  , 

M  <  A/ l!/!l  and  when  Il7ll  <  N  y  ||ij||  <  M 11711  .  Clearly  the 
value  of  s  reflects  exactly  the  properties  stated  above  and 
therefore  is  suitable  for  non-linear  systems. 

2.  The  value  of  s  may  be  used  semi-quantitatively  to 
reflect  on  the  stability  of  round  off  errors.  For  example, 
the  accuracy  of  a  computer  word  is  2"',  L»  in  equation  (2.7), 

A— 2~m  ,  N1  arithmetical  computations  are  required  for 

each  >l  value  calculated.  Therefore,  corresponding  to 

1)711  ,  we  should  have  T» < to  guar¬ 

antee  stability.  Obviously  the  smaller  s  is,  the  more  stable 
is  the  calculation,  and  the  longer  is  the  time  of  stability, 

Tq,  and  the  more  relaxed  are  the  requirements  on  the  computer 
word  length  and  the  step  length. 

3.  For  linear  schemes,  the  generalized  stability  is 

equivalent  to  ordinary  stability  but  they  are  not  equivalent 
for  non-linear  schemes.  Also  if  the  formal  approximate  error 
of  the  scheme  is  O(A'0»  '•  >  K  >  0  ,  then  when  ,  the 

K  th  order  difference  coefficient  of  n  converges  and  is 
bounded . 
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In  recent  years,  Weinberger  et  al.  (see^^-*^ ) 
raised  the  problem  of  the  optimization  of  scheme  and  they 
concentrated  mostly  on  increasing  the  of  linear  schemes. 

In  this  author's  opinion,  generally  a  non-linear  scheme  has 
weak  discontinuities  so  s  is  not  large.  Hence  the  main 

CL 

direction  of  optimization  is  lowering  the  value  of  s.  This 

will  also  improve  the  stability.  Besides  when  s  is  small, 

& 

we  also  have  *  O.  ,  therefore  is  convergent.  For  detail, 
see  [21,  22]. 

We  shall  estimate  the  upper  bound  of  s  below. 
denote  the  computational  error  of  ?»?>,//  .  L,  M,  N,  M 

are  positive  numbers, unrelated  to  h,  t,  a,  b  are  suitable 
positive  numbers,  m»  «  are  positive  numbers  yet  to  be  deter¬ 
mined  . 

m*  -  mix  {2a(2ff  +  X-  l)"',  (20  +  1  +  12JUx,  +  a)(l  +  20)-'} 
m**(«)  -  (80  +  4tf  +  4  -f  48n><,  +  32 •+•  2*1*0)180  +-  4 
+  322  v, a  +  2blv,tf  —  —  *2*  +  I612v,  +  *12*  ]“' 


Condition  (<r,  2,  X,  0,  *)  shows  that  <r»5*(l  —  X)/2  ,  other¬ 

wise  i  <  (80  +  4)»p(16  —  162  -  Ua  +  b  —  Xb  —  2ba)~x  .  Condition 
(a.i,*,0,*)4  shows  that  not  only  the  conditions  above  hold, 

but  also  28  >  m*  when  <x,>(i  —  x)/  2  ,  and  2 $>«**(*)  when 

«r.<(l-X)/2 

The  error  equations  of  scheme  (2.7)  are 

v)  “  i-i(»l •  <?)  ~  X»}  *b  0r>7»,  <p)  —  X» J  +  irtf,,  $>)  —  7il  .  . 

£|(''i,  <f>)  “  L,(f|,  «ji)  —  ?|  '  l 

Detailed  calculation  of  (2.1  +  mr.-),,  I, (17,  «/•))  gives  from 
equation  (1.1)  -  (1.8) 

+  2IHIIJ.  +  -i-0mr*||i}((l|'  -  ±  Xmr'D.O},)  -  2rXA"'(fl,  H'O})) 

1  «  » 

—  —  (v|.,  <jj)  +  5j  S  #)  /  ,  5  \ 

2  /•  1  l»i  \S»£) 
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where 


P,  —  -2B*j(v,  »),  *}.■>,  f,  —  —  mrB*;(»<,  I},,  3,) 

P,  —  ~2ctB**(.v,  i},  q„),  f4  —  —  j}„) 

*.-(f  +  ?)r±U.«i  +  «|) 

X  *  *  '  J 

“  (2*1  +  mr^i,  i(»i  +  £rij,,(£)  +  ?i)  ' 

Pi  —  (2*1 »  X'l  +  fir*},,  9)),  I>»  —  (2i},  JiH  +  fir j},,  9)) 

E,  —  ror(f},,  ;(/}  +  firi},,  9)),  £j  <—  mr(l},,  J(ij  +  Sri],,  9 )) 


We  can  decompose  £,  •»  JS„  +  £«.  .  For  example 

—  2(»} ,  Jii],  $)).  "*  2fir(i},  v(i},,  9)),  £,1-2(*7>y(»l,  9*)),  £,,  —  2Sr(»),  j(ij,,  9)) 

etc.  If  «i  ,  then  we  get  £u  -  ££  +  hi*  from  equation 

(2.4),  etc.  Here 

En  •*  2«,B*j(99/,  i],  l)  *7  2aiB**(f)j<£,  fl,l) 

£*j  2fia,TB**(^9i,  I},,  1)  +  2fia,rB«j(9i*}, ,  f},  1) 

~  2fia,rB*» (*J,  9*},,,  I)  —  23a,r  £«*(*},,  9<}i,  I) 

E?S  -  -2fir(«}„  /(<},  ?)) 

•  •  -  2a,Bit»ifj<pi,  ij,  O  ~  2o»S#;(fl,  9f}„  1) 

It  can  be  proved  that 

|0|r*Hfl.ill‘  <  2 10|  H.ll*  +f  2  )B| 

|22:rfi-'«},  W*Oj))l  <  T  l«K»+  l2*r4^'|f|fyll,  +  W/.Bflll* 

4  )'  ‘ 

.  *  ,  ■  •  V 

U  Xmr'A-  £  iv,ft  <  *Htr*U,V 

1  *  If  I  1  r,  , 

We  shall  first  estimate  the  periodic  problem, 
determine  9  ,  let  9(0*)  •»  coon,  at  some  point  ff'flit,  . 

Theorem  1.  Let  scheme  (2.7)  satisfy  conditions  (»*  4, X, 0, 9) 
and  boundary  conditions  (A).  Then  1.  When  ||7«0*  an<^ 


(3-3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 
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p(.fl,  h,  ht  o.  O  are  not  greater  than  Nh*  ,  then  for  all 

<  T  ^  T#(p)  ,  we  always  have  NCOIIi®.  <  MeLTp.  ,  If  when 

r,A-*o  ,  I!?jII*  —  0(p)  —  oCA1')  then  Tq  is  arbitrary.  Here, 
when  v4  >  0  or  a, «-  a,  ,  *  *  .  If  v„>0  and  also 

a,  —  c,  >  then  /  <  o  •  2.  If  the  conditions  (<r,  X,  X,  6,  0)« 
are  also  satisfied,  and  o,  —  <*.»  ,  then  for  all  h>  p*  , 

the  above  equation  holds. 


Proof:  Equation  (3.2)  holds.  Prom  periodicity, 

-*  .  Multiply  the  second  equa¬ 

tion  of  equation  (3-D  by  V  and  find  the  inner  product, 
from  equation  (1.5)  and  <?(£*)  “  0  ,  it  can  be  proved  that 

there  exists  a  constant  "*•  related  only  to  the  diameter 


R,  such  that 
h 


||$||«  <  m,ll$l!i  <  2mo(ll^ll>  +  \\h\n 


(3-8) 


Prom  the  Cauchy  inequality  and  equation  (3.8),  we  can  get  the 
estimation  equation 


IJf.l  |£,|  <  M ,(||<jll*  +  NIK  +  II7.II*  +  ihV  +  *NIK) 


(3.9)  /136 


If  «,  ,  but  v,  >0  }  then  from  the  Cauchy  inequality, 

lemma  1  and  equation  (3-8),  we  can  get  as  in  [6] 


|2,|  +  |2,|  <  ^  HW  +  e NIK.  +(»  +  —)  W«MI* 

4  '  «*»' 


(3-10) 


|£<l  +  |£«l  < ^ N,ll*  +  elhiiu.  +  w,A-»(i  +  — ) NIl'CNII*  +  117.11')  (3.11) 


If  ,  from  the  periodic  condition  (A)  and  equation 

(3.3),  (3.  *0  we  get  £*  -  2ji  -  2*  -  0;  ;  from  lemma  1  and 

equation  (3*8)  we  also  get 


\iu'  +•  M.NII* 


(3.12) 
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(3.13) 


|-r  II4.U1  +  +  Il7.ll1),  forv,>  0 

I £,*.•]  <  8 

I'f  ||»7.!r  +  A!,fl_,A~J||i;,]*C!l v!I*  +  Ii7.ll1),  Cor  v.  -  0 

Similarly,  we  can  estimate  |£t;,|£,|  •  Substituting 

equation  (3-5)  -  (3.13)  into  equation  (3.2),  we  then  obtain 

+  r  ">011  mil1  +  r(.mO  —  XSt  ,h  —  UXlxt  —  «A/„)||>1.II 

+  IlijlK.  <  ^01,  0)  +  /*(>7.  o)|]'i!ll,  n 


In  this  paper,  let  us  denote 


/*(»!. f)  • -1  +  •-  +  2a  ■+•  (1  —  *ign|<*i  —  a,|)(AM-,||£||}-. 

+  M<.r*!lg<ll7*  +■  sign  v8C II *} II 1  +  Il7.ll1)) 

8(fh  i )  “  (M,jA  '(1  —  sign  v„  +  sign | a.  —  0||)||>}||*  +  Mn*-'(!lif!l7. 

+•  r'tii.liMKIMtl’  +  Ii7.il1)  +  AAsdl-fll1  +  IMII1  +  II7.II1 
+  II7.II1  +  *-111117.  +  r*!|g,ll7J 


Let  *  be  suitably  small,  and  multiply  equation  (3.1*0  by 
x  and  find  the  sum  over  j  «•  o,  |#  2,  ^  —  j  ,  we  then  get 


l4U)Hfe.  ♦  rV.  2  U>).C/)IU.  +  r*  2  (ft  +  2«fl  -  e  ~  JZXAv.  -  «A/, 

-  **  T  2  irOK,).  o)li^(|f.  +  «(.!(/),  0)1 


+  »©m#lO>ll*  +  Cl -0)|| >1(0)11*  +  r  (<r  +  y)  «f)C0)||i, 

+  ripj  ||^(0)||' 


(3-15) 


Let  h  be  sufficiently  small.  If  at>(l  — X)/2  »  then 

take  «•  >  *•  .  Now  coefficients  of  the  2nd  and  3rd  terms  on 

the  LHS  of  equation  (3-15)  are  non-negative  and  may  be  neglec¬ 
ted.  And  finally  we  only  need  to  invoke  lemma  4  and  let 
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iKt)  -  0»Kt)IIU.  *CO  - 


llflCOH!., o 


If  <i,  —  a,  and 
**•  >  C  then  take  «.  ■■  0  .  Otherwise  «,  —  2  .  If 

ff«<0-X)/2  ,  then  take  m  >  .  Because  rIM',  <  8Ai.1||i},||» 

therefore  the  sum  of  the  2nd  and  3rd  terms  on  the  LHS  of 
equation  (3-15)  is  non-negative.  The  proof  below  is  similar 
to  the  above. 

If  further  conditions  (o,i,X,Ot  o),  are  satisfied,  then 
we  may  take  m  —  23  >  m**(o)  .  Prom  condition  (A)  and  equation 

(2.5),  .  After  similar  analysis  we  still  get 

t-i 

(3.1*0,  but  with  Mn  —  M„  ■=■  o  .  The  rest  of  the  proof  may 
follow  theorem  2  in  [6]. 

Note  1.  If  v0 >  o  ,  then  for  the  explicit  /137 

scheme.  If  <*■  —  <»,  ,  then  in  order  to  make  the  principle 

part  of  the  non-linear  error  \t;,  Kt,  ?>))  —  0  ,  and  thus  /«so  . 

In  this  way,  as  long  as  p  does  not  exceed  a  certain  constant, 
the  computation  will  be  stable.  At  the  same  time,  even  if 

some  singularities  occur  in  ,  etc.  on  the  determinable 

01 

long  curve,  as  long  as  we  still  have  the  formal  approximate 
error  Il7.lt  —  «(1)  ,  the  scheme  still  converges.  If  v,  —  o  , 

then  as  long  as  <*,  —  «,  ,  then  for  sufficiently  smooth 

solutions,  t,  2  >  /  —  1  ,  the  scheme  will  also  converge.  This 

explains  why  stability  and  the  2  conservation  properties 
agree  with  one  another.  a,  —  a,  is  optional  for  both. 

Note  2.  For  suitable  ,  the  limit  on  \  may  be 

relaxed.  Suitable  6  may  eliminate  non-linear  error  altogether 
so  that  the  scheme  may  be  used  for  the  global  solution  calcula¬ 
tion  as  well  as  being  perfectly  stable. 

Note  3.  If  the  modified  counter-flow  method  is  used, 
can  be  proved  as  in  [8,  22]. 
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IV.  Two  Step  Method  in  the 

Computation  of  Non-viscous  flow 

Two  problems  exist  for  non-viscous  flow:  1.  When 
0  —  a  <-  o  ,  even  if  a,  ■-  a,,  /,  —  0  ,  we  still  have  Mli  “  .  rMI* 

If  the  gradient  is  large,  then  the  virtual  growth  of  the 
energy  is  fairly  rapid  and  computational  overflow  occurs.  If 
an  artificial  viscosity  term  is  used,  the  resolution  in  front 
of  the  wave  will  be  lowered.  Hence  it  is  desirable  to  design 
a  scheme  that  will  automatically  suppress  this  kind  of  growth. 

2.  In  the  explicit  scheme  of  the  last  section,  but  the 

true  solution  of  the  differential  equation  of  non-viscous  flow 
usually  have  weak  discontinuities.  In  general  /,  <  1  .  There¬ 
fore  we  wish  to  design  schemes  with  •  ^  o  to  guarantee  conver¬ 
gence  . 

Method  I.  Let  n*  be  the  supplementary  value. 

Predictor  v*  “  V  +  tp) 

Corrector  t,  -  i-  (,  +  ,•)  +  <p’>  (4.1) 

If  o,  ~  i  ,  this  is  equivalent  to  the  Maccormack^^ 

method,  but  when  *  «.  a/-  X  ,  9>IU  •  Since 

2 

equation  (4.1)  approximates  equation  (2.1),  (when  p  j  mm  q  ) 
therefore  when  is  increased,  generally  9)11*  also 

increases  and  therefore  the  virtual  growth  of  the  energy  is 
automatically  eliminated,  giving  a  better  result. 

Method  II.  Split  according  to  conservative  type  and 
non-conservative  type.  Equation 

s 

y),  1  —  -i-  C*|  +  O  +  (4.2) 

It  satisfies  the  energy  relation 

M*  -  rlM1  -» f>B* 
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Computational  experience  shows  that  when  J*,J7  is  used, 
the  relationship  between  the  order  of  J*,J~  and  the  direction 
of  wave  propagation  will  markedly  affect  the  computational 
accuracy.  Therefore  for  systems  with  multi-directional  wave 
propagation,  we  should  use  the  average  form,  e.g.  with  J 
substituting  for  J+  etc. 

Method  III. 


“n  +  PrJOi,  <p)  +  T>r Ajj •  <!  —  ij  +  r/G?*,  <p)  (4.3) 

If  j  “j  “  1 ,  r  —  —  hlr~l  ,  this  is  Lax-Wendroff  scheme^-L 

*  4 

If  p—  l,a,  —  l.r-Oi  this  is  Matsuno  sc heme ^ . 

If  p  —  1 ,  a,.—  a,  —  -A-,  r  —  o  ,  then  it  is  conservative  and  gives 

2 

computational  stability.  Sometimes  6  may  be  taken  to  be 
slightly  larger  than  1. 


Theorem  2.  If  in  equation  (4.3),  r  -  0,  a,  —  a,  —  1/2,  p  > 
(3#  +  l)/2,«>0, 


*<  np-  3«  -  +  i)||<p||,]-»#* 

then  when  117,11'  and  0,0  are  not  greater  than  N,  we 

have  #4(01l'<A/ftrp  for  all 

Proof.  Let  the  computational  error  on  the  LHS  be  >u  7< 
Then  we  have  the  following  error  equations: 

A>  *”  /(»!,  y  ■»  v'*  f  7(«j,  #)  +  prJOCij,  <p  +  f"),  q>  +  <p) 

+  Jr . <?),  +  #7  +  PrJUiv,  <p),  f)  +  l 

Multiplying  the  above  equation  with  2i}  i  and  calculating 

the  Inner  product,  we  get 


/138 
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(m  -  l)r|l3,||*  +  7fix\\ 7(3,  <p  +  y)||»-  2  G , 

i 

where 

G,  “  (In  +  m X T;, ,  X*;,  $)  +  pTj(.j(ll,  I p),  <p)  +  ?,) 

G,  -  —  20r(7(>;,  <p),  7(3>  <p  +  $>)) 

G,  —  mr(3,,  7(3t  9  +  $)) 

G,  —  <P  +  <p)«  <P  +  #)) 

C,  •*  ■  rty' »  /(7(*:,  <p),  <p  +  <p)) 

We  can  prove 

1C..  <  M„l  11311*  +  11311*  +  11*11}  +  ||?,||*] 

|G.|  <r*e||7(3,9  +  $)!(*  +  £||$||i 

|C,t  ^  r(3«  +  1)117(3,  «P  +  *)||*  -1-  113,11* 

4(3«  +  1) 

|C,|  <  .r||3.!|*  +  l6fi‘T>a~,i~'tn1\\J(.ij ,ip  +  *)|!*(||$||}  +  ||?>||}) 

•C.i  <  -rll3.il*  +  */„/3*r*#-A-*m*(||vll’  +  MOlflll 

and  hence 

II 3 II ;  r-  r(m  -  1  -  3-  -  i  w*(3«  +  1^“*)I|^II*  +  r(.2p  -  rj  -  3-  -  1 

—  16i*^'«"iwJii<p|l}—  32iV/y*w,m.!l?J!l*)||7(3,  9  +  ^)ll* 

<  3 2 i’w’m.r « " */J * || 3 1| * || 7 ( 3 »  9  +  <p)||*  +  M +  11311* 

+  II?, I!*  -I-  ll?,ll‘  +  r*r*|l3ll«  +  r’4"*llfl|l*ll?jll*) 

Now  let  m  —  6a -i  i,  il7.ll*  be  smaller  than  some  constant,  then 
the  coefficient  of  the  2nd  and  3rd  terms  of  the  RHS  of  the 
above  equation  are  positive.  Finally  we  use  lemma  with 
•"CO  -  I4U)  >* 

Note  1.  We  can  combine  the  weighted  average  conservative 
method  with  the  Pycahob  method,  and  the  Kutler-Lomas-Warming 


method  as  well  as  the  method  in  [13]. 


Note  2.  The  method  in  this  section  is  also  suitable  for 
shock  wave  computation.  The  author  has  applied  it  to  calcu¬ 
late  the  double  shock  problem  of  the  Burgers  equation.  The 
results  show  that  the  two  step  method  can  be  used  for  many 
problems  where  one  step  methods  cause  overflow.  The  calculated 
result  for  «i  ■*  is  the  best.  For  system  with  multi¬ 
directional  wave  propagation,  equation  (4.3)  is  better  than 
(4.1).  The  situation  described  by  Theorem  2  is  far  better  than 
the  Lax-Wendroff  scheme.  For  details,  see  [24]. 

Note  3.  If  equation  (2.9)  is  used  for  non-viscous  flow, 
then  when  r  <  Sh  ,  ,  <  l*» 


V.  Effect  of  Boundary  Value  Error  and 
Boundary  Shape  on  Stability 


Assume  that  on  r,  ,  *P.  f,  represents  the  error 

in  g^.  For  convenience,  let  h  ■*  0. 

Theorem  3.  Assume  in  equation  (2.7)  of  the  boundary  value 
problem  mentioned  above,  the  conditions  Ca,ltX,B,  1)  are 
satisfied.  In  addition,  if  then  >•><),«•  is  bounded, 

therefore  1.  When  ll/ill',  A_,lliMll*  and  h,  h»  Ii*-*,  O  not 

greater  than  AW  ,  for  all  {r  <,  T  ,  we  have 

IlflCOII/#.  <  MeLTp  ,  where  sis  the  same  as  in  Theorem  1;  2. 

If  «i  —  a,  and  the  conditions  C»,  A,  /,<?.! ).  are  also  satis¬ 
fied,  then  when  lliillf,  Nh  the  above  equation  holds. 


Proof.  Multiply  the  second  equation  in  equation  (3.1)  by 
and  find  the  inner  product.  From  equation  (1.5),  we 


get 


fif-llf  1  <  2«;ilW  +  117*11*] 


(5.1) 


i  19 
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where 


rw, 


MP  *** 

vr ||f  (|J 


Therefore  the  results  similar  to  equation  (3.1)  -  (3.13)  all 
hold.  1) 


It  can  be  proved  that 


_  iT'li »¥(,  )  ^  (1  -  e)  2  *Ni-  +  -7)  rAM„||*,||}4 

2  *  **  * 

f,>(i-{)S  *  V  -  +  *«”‘>llf.M« 


*: 


’  P ,  ^  —  ( 2  •'V)  -  —  2  r  «T‘  2  2 

«  *•  "  4  *•  *j  *• 

-  A/,1  r*||f„iu,  +  (*i)'’lli.llh  +  MP  +  MP) 

'V  -  ,r‘  S  -  •  ?  -'i* .. 

2  1  «•  ■*'  2  #j  ■  *;  *; 

-  mmIt*|[|,i|||4  +  r'e-'HMIh  +  Nil1  +  MP1 

f ,  >  2  v'r*  ~  M„(l  +  4a-)t4||fullf. 

2  '  ** 

|2a,B*  *(f]<jPi,  4} ,  I )|  ^  Nil*  ’4"  ^  ’llfttlf*)  . ,  iv 

|2a,B*j(t),  fl,v,  1)|  <  ellflIIJr  +  •  2  *N*  +  «*ICeA»’*>_,ll5.llf* 

*;  • 


(5.2) 


Also  because 

1 5  WH  *  i  £ ■* + -J-  '"iwi’iMif* 

*i  *!  2 

Therefore 


|2u. B»*  (ijf,,  i),  1)||  <  t  2  VY  + 

*  «; 


iiviuiif.iih 


*When  •>.-0  ,  an  additional  term  “•,',lm,li,|}l(Hr+lM,)l  la  added  to 

the  RHS  of  equation  (3.13).  Similarly,  equation  (3.12)  may 
be  estimated. 


|2«,BKfOj,  ffo,  1)|  <  2  V  +  -T 


WasB'iOi,  v>n„  m  <  er'  2  *'d  +  iii.iii.Mi 

*•  *v*h 

|2*«,TB,2(fllt  q-.il,  1)|  <  e  2  *V  +  ^TaillJIfJlyll} 


*; 


ev. 


|2Aa,rB*»(q,  v-^,1,  1)|  «£  2  lP4  +  ilf < i) r«  2  *'*;) 

V'k 

..Aa.rB.jCij,,  *<)„  D!  +  llfwllrjl’lll!. 


/3*40 


Substituting  equation  (3-5)  -  (3.7),  equation  (3.9)  -  (3.13) 
and  the  equations  above  into  equation  (3.2),  we  obtain 

Oi)  +  TmellfJ.H*  +  ximO-UXlVy  -  AAfjX  -  <Ml4)l|iJ,ll‘ 

+■  ll*lll(.  +  f>(fi)*J  2  v’fl*  O  —  5e  —  6«*)  2  ''  i’ 

•:  *: 

f  (7  +  f )  T  (S  *'«*),  <  *<9.  *■>  +  /*0},  j,) II 3 Hi. 


where 


Mi«) 


««• 


—  WwA",||i,||}t(l  —  lignla,  —  o,|) 

Let  e  be  sufficiently  small,  and  multiply  the  above  equation 
with  t  and  also  sum  over  /■■O,  1,  2, •••,{—  1  ,  if  »#>(!  —  *)/2 


9 


then  we  take 
and  if 
If 


™  >  .(» 


Now 


/>.  >0,  f,+  2rn&-0-  12Xly,  -  ,A/„  -  MM  >  0 

is  sufficiently  small,  then  Mil)  >  o 

then  we  take  « >  m**0)  .  Since 


•: 


therefore 

+  (/»,  +  2 me  -  8-  I:*!*,  -  0M,t  -  A/,A;OI|<J,||» 

•  "**  *<£-1  M| i)v  flj  ^  —  A(|iu||f  .  corm 

a  •  * 
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The  rest  of  the  proof  is  similar  to  that  used  in  Theorem  1. 

Note  1.  The  effect  of  boundary  value  error  on  r.  is 
more  pronounced.  If  m  0  ,  then  when  P  <  Nhu,  |j7,ii'  <  AM*  , 

the  explicit  scheme  calculation  is  stable.  If  f^O  ,  then  we 
also  need  to  have  lli.llf,  «£  Nh“+l  .  Also  if  we  take  h,  g,  as 
the  formal  approximate  error,  then  to  guarantee  convergence, 
we  only  require  —  OCA1')  but  on  the  boundary,  we  need  to 

require  llfillr,  *“  0(A*'+l)  .  But  af  =  >h  still  reaches  the 

optimums  as  far  as  error  control  and  the  two  conservation 
properties  are  concerned.  It  makes  s  decrease  by  1,  and  when 
g,  —  o'y/t)  ,  the  scheme  converges. 

Note  2.  If  i  “0  then  when  *’» m  0  ,  the  explicit 

scheme  computation  of  (2.7)  still  converges.  If  g,  o  ,  then 
we  require  J.  ' o, «•  to  be  bounded.  In  the  language  of  numeri¬ 
cal  weather  prediction,  it  means  that  there  should  be  a  smooth¬ 
ing  process  from  the  outside  toward  the  inside,  gradually 
decreasing  in  strength  in  the  neighborhood  of  r*  .  The 
calculation  of  Oliger  and  Sundstroun  proved  this^'*'^  .  /l 4 1 

Note  3.  By  using  methods  similar  to  that  in  [8],  it  can 
be  proved  that  if  a,  —  a,  —  1/2,  v,  >  0  ,  scheme  (2-8)  has  the 

same  stability  property  as  the  implicit  scheme  ( 2.7 ),  There¬ 
fore  it  is  worthy  of  recommendation. 

Note  4.  If  we  calculate  by  using  the  modified  counier- 
flow  method,  it  may  be  proved  similar  to  [8]  that  when 
BAH’ —  o(A‘),  Hi.lifj  —  o(AJ)  ,  the  calculation  is  stable  and  conver¬ 
gent.  But  this  condition  is  too  restrictive.  The  author 
proposes  that  in  the  vicinity  of  T*  ,  equation  (2.9)  should 
be  used  to  be  in  agreement  with  the  transport  properties,  and 
in  the  interior,  equation  (2.7)  for  a,  —  a,  should  be  used  to 
lower  the  s  value.  Shapiro  and  O'Brien,  Williamson  and 


Browning;  had  similar  ideas. 

Suitable  boundary  shape  may  lower  the  estimation  of  the 
upper  bound  of  s. 

Theorem  *4 .  If  the  parameters  in  equation  (2.7)  satisfy 
the  requirement  of  Theorem  3,  R^  is  a  rectangle  and  there  is 
no  error  in  the  boundary  values,  then  even  if  8  —  o  ,  when  ||/1||> 
and  p(<?,  h,  0,  O  do  not  exceed  ,  for  all  <  T  <  T,(p)  , 

we  have  IhlCO'li®.  <  MeLTp  ,  where  if  o,  —  a,  then  *  0.5  , 

if  v,  >  o  ,  then  <<0  }  and  if  “J  and  v#  >  0  , 

then  i  <  — 0.5 

Proof.,  Similar  to  Theorem  3,  but  with  g,  *3  o  .  If 
a*  “  Oj  then  A  —  (fl,  J(ij,  $))  —  0  ,  and  if  a,  ■*?  a,  , 

then  from  lemma  2.3  and  equation  (5.1) 

Ml  +  -1-IH}||||i|ll,ll9ll.llvlli 

2  tv, 

<  -  IHIII.  +  ^  llrjllUOrjll1  +  \\h\n 

2  ev\ 

Further 

$0)1  <  ~  IM.I1*  +  I2m;,«-r||»il|,ll#|ll.tl<pllill«fil» 

When  »» >  t  ,  tIM!I,|]^I1j  <  ,  otherwise  not  greater 

than  A-'A/a/H#}!!*  .  After  the  estimation  by  using  the  above 

method,  we  can  finish  the  proof  as  in  Theorem  3. 

VI.  Effect  of  the  Type  of 
Boundary  Conditions 

Assume  that  on  r  ,  —  g,  where  **>0  .  Its 

bn 
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difference  approximation  Is 


i.(e.  v  +  f  iv(e.  v  +  >w,  p)  -  <•  QiR:  (6.1) 

For  convenience,  let  b  —  X  —  0,  v  ■■  const,  ,  and  the  boundary  shape 

be  sufficiently  regular. 

Theorem  5-  Assume  that  in  scheme  (2.7),  and  (6.1), 
b  —  o,»->o  and  conditions  (<r,  l,  0,  0,  0)  are  satisfied, 
then  when  Mil*  and  h,  in  O  not  greater  than  Nh‘  t 

for  all  <  T  ^  T,  (p)  ,  we  have  lliCOll**'  ^  Meap. 

Proof.  Because 

2fl**(»>,  «,  O  -  2kA  2  “  vh *  2  "■t'* 

•:  *; 

IMIK<con«-(||«||‘  +  «l|«!|{)'> 

Therefore 

l£,l  <  A/«CH4flJ  +  llfill^o  +  ell$!l}.> 

I?, I  <  M1o(r,!|.7,ll*  +  lllollf*  +  et*|M!,  +  f* ll^ii Hr*) 

+  r’lliwll},  +  e||»j||{,  +  l!*illr*) 

IM  <=M»r'(||»7.!|'+  +  «||<i,l|{.) 

Substituting  equation  (3-5),  (3-9)  -  (3.11)  and  the  above 
equations  into  equation  (3.2),  we  then  get 

*  ~m*( '■  .11'  +  r(.m0  —  oM, 4  —  rM*  ~  *A/ ll»l. II1 

+  (2  ~  7f  -  M„e  -  M„e  -  A/*E)||ij||}.  -  ET*||i},||{.(A/. 

+  A/.)  ^  A/H(li i}|J*  4-  I| 5 1|*  +  A-*(||^||*  +  M)l*)ll')ll1 

+  'Ifillf,  +  r'lliJIf*) 

The  rest  of  the  proof  follows  similar  lines  as  Theorem  3. 

^For  proof,  see  another  paper  of  mine.  From  this  get 

|2»,|(*t«1,,)Kc.m|.(H,  +  |*.rf*  +  l-.1V,  ♦•H!-’- 

24 


/142 


Finally  in  lemma  4,  let  —  ||»;(OH4«..  n,  ~  2,  —  o  . 

Theorem  6.  If  in  schemes  (2.7),  and  (6.1),  A>0,»>0 
and  conditions  0,1.  0,0,0)  are  satisfied,  then  when  II/,!1' 
and  p(^,/i,/ a,  ft,  O  not  greater  than  N*4  ,  for  all 

<  7  <  7,0)  we  have 


Proof.  Since  on  Fj 


A.  (A  +  A”)  +  h.  2i?i?,  —  (?*),  —  r(jj,)* 

2  :  .  i 


therefore 


**  **’-  (.1  -  «)  2  (A  +  A')‘  “  M 

1  •: 

p  ^  ^  (fj  +  jj')j  —  ~v -i—  u  •♦* «)  2  (a«  +  ■” w  -hmk 

8  «*  8  *# 

P  ^  Av?!*  y  Cfl  +  ,;■)!  -  m'hi  u  +  *>  2  (A-  +  A*)1  -  wMr,||f„l|f* 

4  £  4 

t,  >  c  i  - «)  S  ca,  t  a;)1  -  Auiiii*.ii/» 

4 


Substituting  equation  (3-5),  (3-9)  -  (3,11)  and  the  above 
equations  into  equation  (3.2),  we  get 

0,.-»(A)  +  rmOHAill'  +  K«0  -  Af  mOIIA'II*  +  20  ”  Oil  All!. 

,  *;‘-o -.)£u  +  o'  +  «=*(.  +  =■) 2 ci  + 

4  *» 

■«  ****  Lm  ~  *  -  2-  -  one  -  «  -  ^  «)  •  2  (A<  +  AO* 

<  AUCIlAll*  MU*  +  *',IIaII’CIIaII*  +  WiO')  +  liiW*  +  r,ll*«ll^> 

The  rest  of  the  proof  is  similar  to  that  for  Theorem  5. 


Note  1.  From  Theorems  5  and  6,  if  &  ■ 0  then  when 
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Il7j!las  p<KV  the  computation  is  stable.  If  i >  0  then 
we  need  to  require  i’ui in  order  to  Guarantee  stability. 

If  we  treat  h,i>  as  formal  approximate  error,  then  when 
||/,||'  =  oU'), ; \i, •-t'U')  ,  schemes  (2.7),  (6.1)  converge.  This 

can  be  achieved  by  suitably  choosing  r*  . 

Note  2.  The  following  externally  imposed  boundary  /l43 

conditions  are  often  used  in  numerical  weather  forecast: 

>?(£>),  c>€/?: 

or 

1<.Q')  ~  WQ)  -  ,(£>••),  qcr; 

where  o** € R, —  KJ  and  distance  h  from  Q.  These  2  conditions 
may  be  considered  as  conditions  (6.1),  with  accuracy  0(1) 
and  0(h)  respectively.  Theorem  5  cannot  guarantee  the  conver¬ 
gence  in  long  term  calculations.  The  calculation  of  Platzman 
and  Mastuno  proved  this.  For  reference,  see  [1*0. 

Note  3*  The  method  in  this  paper  may  be  used  to  prove 
the  suitability  of  many  boundary  treatment  methods  in  aero¬ 
dynamics.  For  example,  #  —  gat  g,  is  often  given  on 

da 

r  according  to  the  famous  Thom  method,  i.e. 

ij(P')  -  S(Q’) - 2 C0C£»  )— +  OU) 

-  -2 t-g,  +  O(A) 

This  is  equivalent  to  the  boundary  value  conditions  of  the 
first  kind  of  n  with  error  o(A)  .  From  Theorem  3,  this 
method  of  treatment  will  lead  to  a  convergent  result. 

Note  4.  We  may  also  adopt  various  boundary  conditions  of 
tp.  For  example,  <p,  —  g,  with  error  g<  .  Let  us  construct 


a  function  5,  such  that  .  Also  let  —  G,  +  <" 

then  —  h  —  AC,  +  i*;,  <?,lr4  —  0  .  Hence  we  can  prove  that  there 

exists  a  constant  m','  only  dependent  on  the  diameter  of 
the  region  such  that 


By  substituting  this  estimation  equation  into  the  earlier 
theorems  as  l!9>|||  ,  we  will  get  the  corresponding  result. 


Note  5-  By  combining  the  various  methods  in  this  paper, 
we  can  prove  the  convergence  of  the  computational  methods  of 
many  practical  problems,  for  instance  the  example  of  a  wind 
tunnel  calculation  on  page  1*40  in  Reference  [1].  Prom 
symmetry,  we  know  <p  —  |  —  o.  on  the  center  line  r,  of  the 
wind  tunnel.  On  the  surface  T,  of  the  object,  we  know 
^  —  o  from  fixed  wall  boundary.  £  is  treated  by  using 
Thom's  condition.  At  r,  upstream,  „  —  „u  .  There¬ 

fore  \p  is  known,  and  {  —  const  by  Daugherty's  method.  On  the 
fixed  walls  of  the  wind  tunnel,  since  v  —  o  ,  ip  is  known. 

£  is  calculated  by  Mueller’s  method.  At  T,  downstream. 


_  ®i_o, 

e»  s. 

This  method  of  treating  boundary  conditions  will  lead  to  a 
stable  and  convergent  numerical  result. 


(,.//////  -  •/>*'//;//////'/' 


7 


.J 


Figure  1 
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VII.  Dynamic  Relaxation  Method  for  Steady  Plow  /144 
Consider  the  problem  of  steady  flow 

OX i  OX}  uTj  OXx  0*\  '  wl/  vXj  '  UXy 

V’4>  -  S  -  /. 

For  convenience,  assume  />  B  0,  v  —  const,  ,  zero  boundary  value 
and  R  a  rectangle.  In  [1,  3,  5]  many  computational  methods 
are  introduced  but  there  is  no  proof  whether  solution  exists 
for  the  corresponding  difference  scheme  or  whether  the  solu¬ 
tion  is  bounded  for  all  h.  In  this  paper  the  following  scheme 
is  given 


(7.1) 


</(»?,  <p)  +  —  —U,  “  ij  (7.2) 

Assume  H  is  the  Hilbert  space  formed  by  the  net  functions 
that  satisfy  the  zero  boundary  value  condition  under  the  Inner 
product 


Iw,  H  -  -jj-  v.J  +  )1  +  “•  S  wv 

The  norm  is  denoted  by  H  •  |!«.  •  Obviously  for  arbitrary 

» €  H,  («',  /i)  .  is  a  linear  functional  in  H.  Hence  there  exists 
FiH  such  that  IF,  w)  —  (»,  /,)  and  1C**',  A)l  <  ||F|j#||n'||/, 

Theorem  7.  In  (7.2),  if  a,  —  a„  >  and  ||F{|„  is  uniformly 
bounded  for  h,  then  for  all  h,  there  exists  at  least  one  solu¬ 
tion  of  equation  (7.2)  uniformly  bounded  for  h. 

Proof.  Multiply  equation  (7.2)  by  w€H  and  find  the 
inner  product.  From  equation  (1.4),  we  get 

—vly,  *> ]  +  C*,  Jiv,  ffO)  ”  —  Oi.  *0 
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But  for  fixed  n  and  <P,  («',  J(*i,  <p))  is  also  a  linear  functional 
in  H.  Hence  there  exists  dytH  such  that  “'J  “* 

(«/,  y(>j,  <p)).  Therefore  the  solution  of  equation  (7.2)  is 
equivalent  to  the  solution  of  the  following  non-linear  opera¬ 
tor  equation: 


ij  -  ±  Un  +  F)  -  q 


(7.3) 


Assume  that  there  is  a  sequence  V*’.  &<pM  —  n{m>  •  When  n~*co  , 
~  9II# -* 0  •  Hence  when  n  is  sufficiently  large,  l|ijulll«, 

is  uniformly  bounded.  Denote  w\  then 

from  equation  (2.4)  and  lemma  2,  we  have 

Take  w  —  Aifm'  »-  Aif'  ,  then  IM(ji<“)  —  V'OIIh  <  M  ,  i.e. 

A  is  a  continuous  operator  and  A/*  is  independent  of  h. 

Now  from  Browder's  fixed  point  theorem,  we  only  need  to  prove 
that  the  equation 


*  —  +  F)  -vO  (7.4) 

is  bounded  for  the  set  of  all  possible  solutions  of  i€lO,l/v) 
Multiply  the  2  sides  of  (7.*0  with  and  find  the  inner 

product  in  H,  we  have 

.  ««  <  All^Ulf  L  <  -J-  U*H,UlF|U 

Theorem  8.  If  2  */~2  <  v\  0,  —  a,  then  the  solution 

of  (7.2)  is  unique. 

Proof.  Assume  n,<p  and  n't?'  are  both  solutions. 

+  +  then 


2  y 


vA»]  +  X«),  <p  +  $)  +  jin,  y)  ■■  o,  Ay  “  fi 

Multiply  the  first  equation  with  n  on  both  sides  and  find  its 
inner  product.  We  also  notice  that 

IhllMM)  <  nvll/.ll*'-.  HyllJ  <  |!Ay||*  <  M‘  <  m;||ij||{, 


therefore 


Hl'lll*  <  2  ||»llli||ij||*||>j||i*  llyll?  llyll* 

<  Zv/Tm^Ihl^lkllJlI^ll.llvll*  <  iy/lv-'MnfM'H  / 

Theorem  9.  Assume  the  ith  equation  in  (7.2)  has  error 
ht  and  the  corresponding  n  has  error  >j  ,  then  under  the 
same  conditions  as  in  Theorem  8,  a/»7£ il/.IJ'  -H  II/jIPJ 

Theorem  10.  Under  the  conditions  of  Theorem  8,  the  follow¬ 
ing  iteration  process  converges  at  least  as  a  geometric  series 
where  t  is  the  relaxation  factor  and  n  the  number  of  iteration: 

,<.♦.> ..  ,<->  +  r[  +  y(,u+",  y(,))  +  /,)  i  ( 7  5 ) 

Ay<*+1>  -  ,<•+'?  +  u  J 

The  two  theorems  above  may  be  proved  similar  to  [8], 

If  —  1  ,  equation  (7.5)  is  then  the  N-S  method  in  [5]. 

Apparantly  this  method  of  calculation  and  proof  is  more  compli¬ 
cated.  In  recent  years,  Davis,  Greenspan  and  Hodgkings  have 
all  tried  to  calculate  equation  (7.1)  with  non-steady  state 
method.  For  reference,  see  [1,  17,  18].  They  call  this  the 
dynamic  relaxation  method.  But  since  equation  (7.1)  is  a  non¬ 
linear  problem,  thus  violating  the  simple  relationship  between 
the  steady  flow  iteration  process  and  the  final  state  of  non¬ 
steady  flow  as  pointed  out  by  Frankel.  Therefore  very  few 
theoretical  results  have  been  derived.  This  author  examined 
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this  problem  from  another  angle,  and  first  proved  that  general¬ 
ized  solutions  exist  for  both  equation  (2.1)  and  (7.1)  under 
very  weak  conditions.  He  also  proved  that  if  the  generalized 
Reynold  number  Re*  —  2v^Tm**i»"*||/,||t,  <  1  then  equation  (7.1)*s 

solution  is  unique,  where 

ItUL  -  Jj,  Vdxxdxu  ||£||L-.  -  Jj, 

lur*.  -ii{iiL:+,ii{iiL.„ 


H*  is  the  closure  of  the  Infinitely  smooth,  finite  subset 
functions  under  11*11* i  ,  and 


mi 


•  mm 


tup  JilL 

IlfllL. 


We  denote  «,  —  2v(»u*)-'(i  -  Re«) 

Theorem  11.  If  £  is  the  solution  of  (7.1),  fW  is  the 
solution  of  (2.1),  both  with  zero  boundary  values,  f  CO  - 15 
{CO— 5  »  then  When  Re*  <  1 

If  COII^C  'TllfCMIL 

Proof: 

§L-  0.1  dX±  ±32  +  +  £?  _  ri>E  „--5L  — 

d»  dx,  0X,  .  .Ok,  ~dft  ..  Qxt  dxt  ..  Qx,  d-r.l  (7.6) 

V»*-*  J 

Multiply  both  sides  of  the  first  equation  by  2 f  and  find  its 
inner  product.  We  get 

•  * 
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From  an  inequality  in  §  1,  Chapter  1  of  Reference  [9], 
ML;  <  IIW-IIL  <  III  Hu  <  «•  lllllu-i.  »  theref  ore 

■J  IIIIIL.  +  — 7  {1  —  Re’Jllflli.  <  o 

ot 


Since  a  systematic  result  has  been  established  for  /lbC 

non-steady  flow  in  our  paper,  we  may  conclude  by  using 
Theorem  11  as  follows: 


2 

1.  For  any  —  «>0  ,  we  may  choose  T»*-—  log 


2||{(0)IU, 


so  that  when  t  >  T.  ,  ||J— f(/)||sS  .  Further  if  we 

assume  that  the  approximate  value  *(#)  of  {(<)  is  calculated 
by  using  any  one  of  the  schemes  in  §  2,  then  we  can  prove 
that  if  conditions  in  Theorem  1-6  are  satisfied  and  <  AM"* 


then  when 


ll«»(T,)  -  f(7»)||  <  e/2 


so  that  ll>/(7,)  —  f||  <  e  ,  i.e.  when  the  above  process  is  con¬ 
sidered  as  an  iteration  process  with  relaxation  factor  r,  the 
process  is  convergent. 

2.  Iterative  convergence  depends  on  many  factors.  If  v 
Is  too  small  or  ll/Ju  too  big,  the  conditions  of  Theorem  11 
will  both  be  violated  so  that  the  total  iterative  procedure 
may  not  converge.  If  h  is  too  big,  then  no  matter  how  big  is 
Tq,  j(0  will  not  converge  In  {(<)  .  If  t  is  too  small, 

the  amount  of  work  will  be  increased;  if  t  it  too  big,  then 
the  conditions  of  Theorems  1-6  will  be  violated  and  the 
iterative  process  may  not  converge.  From  the  process  of  proof 
we  see  that  the  iterative  convergence  is  related  to  initial 
iterative  error,  boundary  shape,  type  of  boundary  condition  as 
well  as  to  boundary  value  error.  This  agrees  with  the  exper¬ 
ience  described  by  Greenspan^^ ,  Roach^^,  etc. 
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3.  The  amount  of  work  may  also  be  estimated.  For  example 
if  equation  (2.7)  is  used  in  the  calculation,  the  total  number 
of  iterations  is  about 


2||!(o)IL 


4A7  \i,  V* 


It. 
I  "•! 


2|IIC0)k 


Apparently  when  v> 8  are  large,  Kq  is  small.  When  U/»||,  |||(0)||t, 
or  M,  {^(gradient  correspondong  to  the  solution)  are  large, 

Kq  is  large. 


When  Theorem  11  is  applied  to  local  district  numerical 
weather  forecast,  it  shows  that  if  v  is  fairly  large  and  the 
n  values  on  the  boundary  are  fixed,  then  when  A -*0,r -*■<>, 

*(0  ~*fO)  and  when  TQ  is  sufficiently  large,  ff  (0 
approaches  £ ,  i.e.  the  forecast  value  rapidly  transits  toward 
equilibrium.  Therefore  the  sponge  boundary  conditions  in  [20] 
should  be  used.  For  more  detailed  proof  and  application  of 
some  of  the  results  in  this  paper,  see  Reference  [21-2*1]. 
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Abstract 

There  are  five  important  subjects  in  the  theory  of  difference  methods:  1.  eon- 
struct  ion  of  scheme;  2.  Stability  of  scheme;  3.  Computation  of  flow  without  viscous 
term;  4.  Interference  of  boundary  conditions;  5.  Solution  of  steady  flow.  A  lot  of 
work  has  been  done  in  this  field**-*1,  but  no  systematic  theory  has  been  developed. 

In  this  paper,  systematic  results  are  provided  with  particular  reference  to  two- 
dimensional  vorticity  equation.  Two  classes  of  schemes  are  constructed  based  on  con¬ 
servation  and  transport  property.  A  new  concept  concerning  generalized  stability  and 
optimization  of  nonlinear  scheme  is  introduced.  The  index  s  to  generalized  stability 
for  periodic  problem  is  estimated,  which  shows  the  relationship  of  stability  and  con¬ 
servation.  Two-level  schemes  are  given  to  decrease  the  index  a.  The  influences  of 
bound  shape,  boundary  conditions  and  their  errors  are  fully  discussed.  Finally  the 
existence  of  the  solution  of  steady  flow  is  proved  and  a  dynamic  reluxation  method  is 
proposed.  Several  examples  are  given,  showing  that  the  methods  suggested  here  are 
useful  both  for  numerical  weather  prediction!  and  in  the  field  of  aerodynamics  ..uJ 
other  subjects. 
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ON  THE  KERNEL  FUNCTION  COLLOCATION  METHOD  IN  STEADY  SUBSONIC 
FLOW  FOR  WIND  WITH  CONTROL  SURFACES* 


Chen  Jing  Song 

Nanjing  Aeronautic  Institute 
ABSTRACT 


In  this  paper  the  forms  of  the  lift  distribution  function 
for  wings  with  control  surfaces  in  steady  subsonic  flow  are  anal¬ 
yzed.  Methods  for  treating  kernel  singularities  of  linear  inte¬ 
gral  equations  and  singularities  of  the  lift  distribution  function 
are  discussed  and  the  numerical  solution  to  the  integral  equation 
is  given.  The  method  in  this  paper  may  be  used  to  calculate  the 
lift  distribution  of  the  wing  by  itself  or  of  the  wing  with  par¬ 
tial  or  full  control  surfaces  on  both  the  leading  and  trailing 
edges.  The'  numerical  values  agree  well  with  experimental  data 
and  have  the  same  accuracy  as  those  given  by  other  theories.  As 
compared  with  the  vortex-lattice  method,  the  computational  storage 
requirement  for  our  method  is  smaller  and  the  computational  time 
is  less.  The  numerical  computation  may  be  carried  on  smaller  com¬ 
puters  . 

I .  FOREWORD 

At  present,  the  numerical  methods  to  solve  the  lift  surface 
linear  integral  equation  are  categorized  into  two  kinds:  The  vor¬ 
tex  lattice  method  and  the  kernel  function  method.  The  advantage  of 
the  vortex  lattice  method  is  that  it  is  not  necessary  to  assume  a 
priori  the  form  of  the  lift  distribution  function,  and  hence,  in 
principle,  there  does  not  exist  any  difficulty  to  compute  problems 
Involving  partial  control  surfaces.  However,  to  achieve  a  certain 
computational  accuracy,  the  number  of  vortex  lattices  should  be 
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large  and  this  will  certainly  result  in  a  higher  order  for  the 
linear  algebraic  equations  of  the  coefficients.  The  advantage 
of  the  kernel  function  method  is  that  it  requires  very  few  coll¬ 
ocation  points  and,  therefore,  the  order  of  the  equations  to 
solve  for  the  undetermined  coefficients  is  low.  Its  difficulty 
lies  in  the  fact  that  the  lift  distribution  function  with  match¬ 
ing  boundary  conditions  and  edge  conditions  must  be  chosen  a 
priori.  Reference  [3]  made  use  of  a  method  of  approximate  expan¬ 
sion  to  study  the  pressure  amplitude  characteristics  of  the  lead¬ 
ing  and  trailing  edges  as  well  as  the  corner  when  the  control  sur¬ 
faces  oscillate  harmonically  and  obtained  the  "form  for  the  control 
surface  lift  distribution",  thus  making  it  possible  to  solve  the 
control  surface  problem  with  the  kernel  function  method. 

There  are  two  methods  of  treatment  when  the  kernel  function 
method  is  used  to  solve  the  control  surface  problem.  One  of  these 
is  to  substitute  directly  into  the  integral  equation  the  control 
surface  lift  distribution  as  obtained  by  the  approximate  expansion 
method  to  find  the  downstream  velocity  induced  at  the  collocation 
points  on  the  wing  (including  the  control  surfaces).  The  equiva¬ 
lent  downstream  velocity  Is  then  obtained  by  subtracting  this  in¬ 
duced  velocity  from  the  actual  downstream  velocity  at  this  colloca¬ 
tion  point.  In  this  way,  we  may  solve  for  the  lift  distribution 
that  satisfies  this  equivalent  downstream  velocity  with~the  kernel 
function  method,  while  regarding  the  wing  with  control  surface  as 
one  without.  Then  by  super-imposing  this  distribution  to  the  con¬ 
trol  surface  lift  distribution  mentioned  above,  the  lift  distribu¬ 
tion  on  a  wing  with  harmonically  oscillating  control  surfaces  is 
obtained.  This  kind  of  method  is  called  the  equivalent  downstream 
velocity  kernel  function  collocation  method .  ^ 

To  obtain  a  stable  solution  for  the  equivalent  downstream  velo¬ 
city  kernel  function  collocation  method,  the  number  of  collocation 
points  should  be  somewhat  large.  It  is  also  necessary  to  select  an 
expansion  factor  for  the  control  surface  lift  distribution,  but  it 
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is  difficult  to  make  it  satisfy  the  edge  condition  of  a  flat 
plane  while  keeping  the  second  order  derivative  of  the  lateral 
load  continuous.  Shortly  after  the  completion  of  the  work  in 
this  paper  in  August  of  1 9 7 3 »  we  saw  that  some  progress  had  been 
done  along  this  direction  in  [6].  Some  of  the  problems  are  solved 
but  the  region  of  integration  has  to  be  subdivided  into  many 
smaller  regions  during  numerical  integration.  The  technique  of 
using  just  a  few  collocation  points  and  integration  points  as  in 
the  control  surface  kernel  function  method  cannot  be  used.  In 
this  paper  we  use  the  control  surface  kernel  function  method  to 
determine  the  lift  distribution  of  the  control  surfaces.  Limited 
by  space,  only  the  results  for  the  steady  flow  will  be  presented 
in  this  paper. 

II.  FUNDAMENTL  EQUATIONS 

The  fundamental  equation  of  steady  subsonic  lift  surface 
theory  is 

^  »•»  < ! ) 

where  w  is  the  downstream  distribution  function  over  the  wing; 

p  and  are  the  incoming  flow  density  and  velocity;  AP  is  the 

lift  distribution  function  of  the  wing  and  is  a  constant; 

K  _ _ 1  ■  [|  -p  *  ~  |  1  is  the  kernel  function;  Mm  is 

(y-1?)'  I  \/(*  -  s y  +  P‘(.y  -  I?)1  J 


the  incoming  flow  Mach  number,  p  —  vl—  Mi. 

According  to  Figure  1,  by  transforming  the  coordinates  in 
Equation  (1)  and  making  use  of  the  conditions  that  AP  is  zero  at 
both  tips  of  the  wing,  we  obtain 


where  w  -  AnpU.wii,  y),  *  -  I)’  +  pKy  -  f)Y, 


.  In  the  above,  the 


dimensionless  coordinates  are  respectively 
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t  „  „  (*\  (JL\.  (±\  ( 
r*  >  *  £*  1,  )*  l,  J»  l,  )>  l,  ) 

(i>  hm-  © 


physical  plane 


transformation  plane 


- 

2'h 
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t 

Figure  1.  The  physical  plane  and  transformation  plane 
of  the  wing. 

The  meanings  of  the  other  symbols  are  shown  in  Figure  1, 


To  obtain  a  convergent  solution  with  suitable  accuracy,  to 
reduce  the  required  machine  memory  and  to  shorten  the  computational 
time,  before  solving  the  integral  Equation  (2),  we  must  first 
select  a  suitable  lift  distribution  function,  find  the  best 


positions  of  the  integration  point,  the  downstream  point  and 
use  a  stable  and  effective  numerical  integration  technique. 

III.  SELECTION  OF  THE  FORM  OF  LIFT  DISTRIBUTION  FUNCTION 

According  to  linear  theory,  the  lift  distribution  may  be 
divided  into  two  parts:  the  normal  part  AP^  as  produced  by  the 
incident  angle  when  your  angle  is  zero  and  the  subsidiary  part 
APcsr  as  produced  by  the  zero  angle  when  the  incident  angle  is  0. 

AP(|,  5)  -  3)  +  5)  (3) 

The  function  APR  is  required  to  be  continuous  on  the  wing 
with  an  inverse  square  root  singularity  at  the  leading  edge  of 
the  wing,  to  satisfy  the  Kutta-Joukowski  condition  at  the  trailing 
edge  and  to  vanish  at  the  side  edges  with  a  square  root  form. 

Hence,  the  tangent  direction  of  APR  may  be  chosen  to  be  the  first 
few  terms  of  the  accurate  thin  wing  solution  and  the  lateral 
direction  to  be  the  elliptic  distribution  solution  of  a  (long, 
narrow  slender  body)  wing  [7,8]. 

^*<5.  a>  -  Wi  -  v1  u.c3> y^-rj  +  *‘(3)<i  -  vr*  cl<) 

where  V*(3)U  -  ey*g  +  ••;] 

f.(3)  —  +  a.iU i(»j)  +  amiUt (j)  +  •  •  .(n  —  0,  1,*  •  •,  N~)  (.5  ) 

(N  +  1)  is  the  number  of  terms  along  the  tangent  direction  of  the 
lift  distribution  function;  Um(n)  is  the  Chebyshev  polynomial  of 
second  kind;(m  =  0,  1,  ...,  M)  is  the  number  of  terms  in  the  later¬ 
al  direction;  anm  are  the  undetermined  coefficients  to  be  found 
by  satisfying  the  boundary  conditions.  It  is  easy  to  see  from 
Equation  (4)  that  AP^  does  have  the  proper  edge  properties. 

The  zero  incident  angle  slender  body  lift  distribution  [9]  of 
the  partial  control  surface  (with  axis  of  rotation  along  the  lead¬ 
ing  edge)  is 
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where  Ap  is  called  the  control  surface  lift  distribut ion  func- 

cs 

tion.  It  is  proportional  to  the  logarithmic  term.  Ap^  is  called 
the  normal  lift  distribution  function.  Similarly,  the  zero  inci¬ 
dent  angle  lift  distribution  function  of  the  wing  with  partial 
control  surface  consists  also  of  two  parts, 
similar  to  Equation  ( h ) 

form  supplied  by  [3]  in  conjunction  with  the  logarithmic  term 
solution  of  the  wing  type  discussed  above  as  [73 


The  form  of  AP^  is 

AP  may  be  expressed  according  to  the 
c  s 


where 


APciCf,  3)  *“  ~^-y  Si>V  1  —  31  [£<.(3)  +  £1(3X1  +  J)  +  •  •  •  ]<pcj(£,  3) 


,  {[I  —  tfc.  +  0  -  Wl  -  +  E -  E, 

9>cX§,  3)  ■=  2~,!n  - 

■  -I  1 


(6a) 


(6b) 


Kf  -  ic.y  +  £iJ'/J  -  E. 

E  -  —  3)>  3  ^  0 

. >— ^(3  +  3i/)(3».  P  3),  3  <  0 

Summarizing  the  above  discussions,  the  lift  distribution  function 
of  the  wing  with  partial  control  surface  at  zero  angle  is 

APc«(|,  3)  “  ~^yy  S0y/\  —  3*  { f 2 >  +  £1(3X1  +  £)  +  •••l<po(f,  3)  ^  j 

A  fv/j 

+  l£yc(3)(^j-^;-|J  +£i/t+u(3)(l  —  |,),/*  +  £i/c+j)(3)(l  —  +  •  •  •  1} 


In  Equations  C6a)  to  (7),  ICS  is  the  number  of  control  surfaces 

that  move  together  in  the  right  half  of 
the  wing;  for  the  meaning  of  fc.,  31. 
and  3»*  ,  see  Figure  2;  JC  is  the  num¬ 

ber  of  tangential  terms  of  AP  .  For  the 

c  s 

case  of  leading  and  trailing  edge  control 
surfaces  that  rotate  together,  JC=2. 
Otherwise  JC=1.  For  the  tangential  terms 
of  APr,,  at  least  two  terms  should  be 


Figure  2.  Diagram  show¬ 
ing  the  control  surfaces 
on  the  wing 
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taken.  It  Is  not  difficult,  to  verify  from  Equation  (6)  that  AT 

0  s 

has  the  property  as  pointed  out  in  [3]:  namely,  that  it  is  zero 
at  the  leading  and  trailing  edges  of  the  wine,  finite  at  the  side 
edge  of  the  control  surface  and  continuous  on  the  wing  surface 
outside  the  control  surfaces. 


IV.  THE  BEST  POSITIONS  OF  DOWNSTREAM  POINT  AND  INTEGRATION 
POINT  AND  TECHNIQUE  OF  NUMERICAL  INTEGRATION 


It  can  be  proved  that  to  minimize  the  difference  between  the 
approximate  values  of  the  low  speed  thin  wing  lift  and  the  total 
lift  and  their  respective  exact  values,  the  tangential  downstream 
point  and  the  lateral  downstream  chord  should  be  respectively  [10] 

*' "  “*  (2F+1)  V-’.J.-.JO 


(8) 

(.9) 


where  I  is  the  number  of  downstream  points  of  each  chord;  R'  is  the 
number  of  downstream  chords  of  the  right  half  wing.  The  best 
ratio  of  R'/I  may  be  obtained  following  [7]  for  various  wing 
shapes  and  M^. 


From  the  functional  forms  as  shown  In  Equations  (4)  and  C7 ) , 
we  can  use  the  Gauss-Mehler  quadrature  with  \/(l  —  jp/(l  +  f) 
as  weight  to  maximize  the  accuracy  of  the  tangential  numerical 
integration  of  Equation  (2) 

«*«*>-  UOa) 

where  J  is  the  number  of  integration  points  of  each  chord.  The 
integration  point  and  weighting  coefficient  are  respectively 


H, 


-co  .[« 
127 


-  1 


■h  1 


7* 

2J  +  1 


O  -  {,) 


ClOb) 
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If  f(£)  is  a  polynomial  with  order  less  than  (2J  -  1),  then  Equa¬ 
tion  (10a)  is  the  exact  value. 


Substituting  Equation  (4)  (or  Equation  (7))  into  Equation 
(2),  we  may  obtain  the  general  form  for  the  lateral  integral 


“  J-.  v/rr3rG(i;’  d3 


G(y_,n)  is  the  tangential  integral  of  Equation  (2).  We  may  assume 
that  it  is  an  analytic  function  of  n.  Since  it  is  impossible  to 
numerically  integrate  Equation  (11)  directly,  we  first  expand 
G(y_,n)  into  a  Taylor  series  and  then  integrate  termwise  to  get 

«,)  --±n.  e> 

/  - 1  1  sjf  Jil 

+  }  (.12a) 

where  S  is  the  number  of  integration  chords  of  the  whole  wing; 
and  Hs  are  respectively  the  integration  point  and  weighting  coeffi¬ 
cient  of  the  Gauss-Mehler  quadrature  Equation  (11)  with  l/\/  1  —  j1 
as  weight. 


3-“  “  O-  1,2,  •••,5)1 

H,  -  n/S  I 


(12b ) 


It  will  be  proved  below  that  Equation  (12a)  may  be  greatly  simpli¬ 
fied  by  appropriately  collocating  the  positions  of  the  downstream 
chords  and  the  integration  chord. 


Since  the  positions  of  the  integration  chords  mentioned  above 
•••>  are  zero  points  of  Tg(n),  therefore 

T.Cj)  -  2J-  n<2  -  3<>  QL3a) 

The  logarithmic  derivative  of  (13a)  is 

n(3)/T,  <3)-S-z—  a3b) 

3  3* 

If  we  let  2R ’  +  1  *  S,  i.e.,  take  as  the  zero  point  of  Tg(n) 
then  n  and  v  become  crossed.  Also  from  Equation  (13b) 

”S  * 
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V  _ ! _ 

X'  ~  V‘ 


(13c) 


Hence,  Equation  (12a)  may  be  simDlified  as 

k>.)  ~  -  j  £  s-  . 

.  (13d) 

+  \  S ; - ^7] 

1  ^  ,-j  (*,  —  2,)'  1 

We  shall  prove  below  that  the  square  bracket  portion  of  Equation 
(13d)  is  equal  to  S.  Change  the  variable  n  in  Equation  (13b)  to 
y  and  differentiate.  We  then  get 


S  7 - 3— s-  “  r»’(z)/T,(z) 


C13e) 


According  to  the  property  of  Chebyshev  polynomial, 

O  -  iWX*)  -  ,n(z)  +  -  0 


Cl  —  jf^rtCj;)  —  +  S'r,^)  —  o  (13f ) 

Therefore,  we  get  from  Equations  (13e,f)  and  the  fact  that  y^  is 


the  zero  point  of  T'(y): 


v  — ! - 1! — 

a-tf) 


(13g) 


Since 


(1  —  3?)  —  (1  —  *!)  +  —  3,)  —  (*,  —  3,)',  therefore  , 

±  ,-l^rr  -  (1  -  ±T~i—  +  2t,  2—1 - ±, 

1  Cjfr  JJ<)  «■*  (jff  —  Ji)  #»l  Jff  —  J<  I  •  I 


(13h) 


It  may  finally  be  proved  from  Equations  Cl3c,g,h)  that  the  square 
bracket  portion  of  Equation  (13d)  is  S.  Hence,  the  numerical  inte¬ 
gral  (12a)  is  finally  simplified  [10]  as 


*  v-«  O  ~  5<')G(?,,2,)  . 

“  -7  £  + -sc<'- 1-' 


If  G(y,n.)  is  the  (2S  -  1)  or  lower  order  polynomial  of  n,  then 
Equation  (1^)  is  the  exact  value. 


V.  INTEGRATION  OF  THE  INTEGRAL  EQUATION  AND  ITS  SOLUTION 


1 .  Case  of  grille  a  X  0  and  yaw  angle  6  =  0. 

Substituting  the  lift  distribution  function  Equation  (4) 
into  the  integration  formula  (2)  and  applying  Equation  (1*0,  we 
get 

!?•!',  r-l  --  »  ± »  - +  t.) 

VUL  s  (*,  -  2<)  -  (15, 

where 


f1  r  x  —  f 1  ( ) 

CX*.  2>  “  ],,  *X£,  3>  1 1  +  — jp] 

2)  “  (£..(3)  +  £.(3X1  +  f)  +  +  f)  +  •'•1Vr+' |  (_17) 


In  Equation  (16),  the  function  being  Integrated  has  a  jump  discon¬ 
tinuity  at  ?  =  x  along  n  =  y.  Because  of  this.  Equation  (16)  is 


rewritten  as 


CX*.  3)  -  1<*>X|,  3)  -  <AX|,  *)1  [1  + 


Thus,  the  singularity  of  the  Integrand  in  the  first  integral  may 
be  found  by  using  the  quadrature  Equation  (10).  Jump  discontinuity 
still  exists  In  the  second  Integrand  which  is  difficult  to  inte¬ 
grate  analytically.  Hence,  we  expand  the  analytic  function  <#>*(§, j») 
into  a  Taylor  series: 

4>Xf ,  2)  —  $X*.  *)  +  ({  —  +  •  •  •  (19 ) 

while  from  Equation  (2a)  and  Figure  1,  it  can  be  shown  that 

<*-»>- l  VK*)K«  -  ff)  +  l*/Ki)H«  -  MOl -«  +*>  (20) 

Substituting  Equations  (19)  and  (20)  Into  Equation  (18),  the  second 
integral  can  then  be  integrated  analytically.  If  only  the  first 
two  terms  of  the  series  are  taken,  then  after  some  algebra,  the 
result  of  integration  is 
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where 


•  1  <•! 

•  jl  +  X~L~^'  1  +  2  *)</'»,(*'»  y»  «l)  (21) 

1  A  tt  J  *-0 

y)  *“  lfr(y)  +  f.(y)(l  +  r)  -f  f>(y)j(l  +  t)  +  •••]  (22) 

4>#i(*,  y)  —  —  l+mCs,  y)]  (23) 

*»(*,  y.sO  -  {'U*  ,y,  l>“jrj  2  O  -  JJ)w[l  +  -1-]}  (.2*0 

^.(i,  y. «?)  -  -  O  +  *)]■/'  «(*,  y, »?) 

+  -  tv -  iff\  &  <' -  era  -  « [•  ^11  (25  ’ 

y,  fl)  -  2  +  T7*^(KJ— ii,)*  +  /*’(y  -  +  Wy-^T1)  (26) 

‘><*1 J 

/t, a,  y,  o--j  {(s  - 1*)1  -  a. -  ?y  +  a  -  6*.)i<5  -  £*.v 

+  ^*Cy  ~  O']'"  +  (£,.  -  *)[(£..  ~  *)*  +  A'(y  -  ?)*]"»  -  /?'(y  -  *})‘ 

1 1(6*.  -  O'  +  PK9  -  0‘Y/‘  -  a,  -r)i!  Ldl  1 


CXy,  y)  “  limG*(y,  tj) 
3-1 


Substituting  Equations  (21)  and  (28)  into  Equation  (15) »  the  down¬ 
stream  distribution  expression  for  a  /  0 ,  S  =  0  is 

E&i, it)  _  .  .-IzL-  y  y  ILull 

s(2y  +  l)  ;r|  Cy.-a*)' 

U.(a*Xi  —  fi)  +  fiCa-Xi  —  $!)  +  •••) 

,  [i  +  *<  ~~  >i  1  t  *»*£.  C29) 

1  Vi*  -  t,y  +  m.-  n.y]  v  +  'ft 
lf«<y,)0  —  |  <)  +  fiCy.Xi  —  {<)  +  •••) 

■  |l  +  p^M  -  f  ±  }--%  f.  j.) 

*  |jr*-{,N  s  pi  (y.-a *)V 

+  *svM(.ti,  y.,  y,) 


M5 


Tho  last  two  terms  in  the  above  equation  are  the  correction  terms 
that  appear  when  the  jump  discontinuity  in  the  integrand  is  elim¬ 
inated.  In  the  above, 

I 

VcHii,  y,  3)  -  2  y) '  *»'(•*>  y>  O  ( 30 ) 

r-o 

The  case  of  a  =  0,  3  ^  0. 


From  the  simultaneous  Equations  (2), 

ItT'X*.  V.)  _  _  H  V* Ct  ~  q])GCsk(.X'>  3<) 
ApUi  "*  S  Tr[  C y,  -  3.)' 


(7)  and  (1*0, 

+  *SGcsn(Xrf  }fr) 


where 


GciXy,  3)  "•  Gcsn(x>  3)  +  GcshCx,  3) 


we  obtain 

C3D 

(32) 


Ccim  “  j_(  UX3)  +  ^1(3X1  +  f)  +  •••]<p«(S>3)[l  +  --  (  33a) 

Gan  —  (£7X3)  +  f<;c+i>(3)(l  +  |)  +  +  ~  ^  (3**a) 

Equation  (34a)  is  similar  to  Equation  (16).  The  result  of  inte¬ 
gration  is 

,  J 

GcsniXt  3)  *"  t  *  ,  2  ({7X3)0  —  {>)  +  {(70+0(3X1  ~  |y)  +  "'l 
•7  +  I  f«| 

•  [1  +  H  +  2)  1 i)^*/**,  y,  <j)  C3*ih) 

Except  for  the  jump  discontinuity  at  n  =  y,  (  =  x,  the  integrand 
in  Equation  (33a)  has  a  logarithmic  discontinuity  at  the  leading 
edge  of  the  control  surface.  The  method  of  eliminating  the  former 
is  as  described  in  the  above  section.  The  logarithmic  singularity 
is  eliminated  by  a  method  similar  to  that  used  to  eliminate  the 
Jump  discontinuity.  To  make  the  correction  term  that  changes  the 
integral  from  singular  to  non-singular  consistent  along  the  lateral 
direction,  the  lateral  part  of  the  non-singular  term  is  also  cor¬ 
rected,  The  final  result  of  Equation  (33a)  after  this  treatment  is 
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2*  --A 

Cc^(i,  S  I f.CgXl  -  {,)  +  *,(3)(1  -{))  +  •••] 

’  Vi  )  “^fx  !/)| 1  +  /^“'J  +  VcciCf/,  a,]) 


* 

^to(x,  2,  ’j)  —  v  [*/>cs<,(jr>  *)</>,,(*,  y,  3) 

ICS 

~  <t>ic.cs,(.(c.,  q)4>cs.,(y,  »j)  I 

•  ■  1 

4cte(i,  a)  -  ftf/y)  +  f,(y)(l  +  s)  +  •••]<*>«(*,  y) 

*«'(*»  y)  -  d^cs„(St  ^/Qg 

4>ci*>(.y,  3 )  -  {L0(fc.,  3)  -  (i  _  |.)w 

•blKSi  ~sec),  + 

4’cs.Xy,  fl)  -  {L,(|c.,  3)  -  g  (1  _  £)'"(!,  _  {c>) 

*  «“l +  £i)1/J  -£.|} 


(33b) 


(35) 

(36) 
C37) 

(38) 

(.39 ) 


3)  -  (1  -  |c.)ln|  1(1  ~  |c.)»  +  £«]./>  _  £j 
+  O  +  fc.) ln|  [(I  +  fc.)'  +  £*]•'»  _  E.|  _  £ * 

• |n  jilli  *  fr*)1  +  +  O  +  gc.DUO  —  fc.)»  +  Ej]1^  +  (1  -  ^c>)j 


2  (40) 


L,(Sc-  y)  ”  S<=-  +j{0  -  {cJ'Inl (Cl  -  fCt)>  -I-  Ei]1'1  -  E.l 

“(1  +  |c.)ln|f(l  +  §c.)‘  +  £•>)*'»_  £,|) 

~  —E.ilO  -  $c.y  +  —  [(1  +  £c,y  +  e’]'*} 

+*ct.(fc,  3)  -  +  f,<3xi  +  £)  +  •"][]  + 


3)  -  [—•*««($,  5)] 


H-fc. 


(41) 

(42) 

C43) 
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The  expression  for  downtream  distribution  of  a  =  0,6  /  0  is 
obtained  by  using  Equations  (3D,  (32),  (33b)  and  (34b): 


Wcnixi.  D  „  _  m  -  y  y  1  ~ 
ApV'.  SdJ  +  O  frj  fzi  (*,  -  2,)1 

•  {UoC3,)0  -  p  +  f.ca-xi  -§))  +  •••) 


•  ^  <pcs(.§,,  2-)  +  Igjcij'Xl  —  §,)  +  ^yc+i>(2«)0  ”"{})+  •  ’  ‘  J) 

.  Ii  +  u  1»!I-  v* 

l  ’Sa,  -  iy  +  PCy,  -  P,J  2^  +  D-. 

UftCjfrXl  —  £/)  +  SlC^rXl  —  fj)  4-  _  £() 

•  TciCfn  l>)  +  UvcCpO  {,)  +  fuc+uCpO  —  fi)  +  *  •  *  1) 

•  [i  +  T--~y  -  ~  2  51  ~  t„  2-)  +  fcXii,  ir,  2'X 

+  *$[«foc*C**»i>,i>)  +  i>,  *,)] 


The  un-explained  symbols  in  this  section,  with  the  exception 
of  lgjcW+guc+i)(xKl+x')  +  ---]  ,  are  identi¬ 


cal  to  those  in  the  last  section. 


3.  Determination  of  the  undetermined  coefficients  a»* 


We  write  Equations  (.29)  or  (.*>*0  as  a  matrix  equation: 

~  »  (45) 

where  {«.„}  and  are  column  matrices;  is  the  boundary 

condition  at  the  optimal  downstream  point  (x^,  iLp)»  i.e,,  the 
inclination  of  the  object  surface;  [(A.),!  is  the  matrix.  The 
elements  of  the  matrix  may  be  obtained  from  Equations  (29)  or  CM), 


After  anm  are  solved  from  Equation  (45),  the  lift  distributions 
and  the  other  aerodynamic  quantities  of  the  wing  surface  and  the 
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control  surfaces  under  given  conditions  may  then  be  found  from 
Equations  (4)  or  (7). 

VI.  SAMPLE  CALCULATION  AND  CONCLUSIONS 


The  above  computational  method  has  been  coded  into  a  program 
for  the  X-2  computer.  Calculation  shows  that  the  numerical  re¬ 


sults  of  our  method  agree  well  with  experimental  data  with  the 


same  accuracy  as  other  theoretical  values.  Figure  3  shows  one 


sample  calculation,  which  requires  only  two  minutes  of  machine 
time  on  a  709  computer. 


2 

atei+tt 

♦  [2]  =  8,  /  “  10 
—  [7]  K  “12,  I  “1 
12.  3 


Figure  3-  Comparison  of  the  computed  values  in  this 
paper  to  the  experimental  data  and  other  theoretical 
values  <*/„ •.<>.«,  3-<m«) 

Key:  1 — experimental  data;  2 — theoretical  calculation; 

3 — ratio  of  lateral  to  tangential  dimensions; 

A — trapezoidal  ratio 


From  the  limited  numerical  results,  it  has  been  demonstrated 
that  the  method  presented  in  this  paper  is  reliable.  Compared  to 
the  vortex  lattice  method,  our  method  requires  only  a  few  down¬ 
stream  points,  thus  greatly  reducing  the  machine  stroage  require¬ 
ment  and  computational  time.  The  numerical  computation  may  be 
carried  out  on  small  computers. 
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ON  THE  KERNEL  FUNCTION  COLLOCATION  METHOD 
IN  STEADY  SUBSONIC  FLOW  FOR  WING 
WITH  CONTROL  SURFACES 

Chen  Jing-song 

(Nanjing  Aeronautical  Institute) 

Abstract 

In  this  paper,  the  forms  of  lift  distribution  function  for  wings  with  control  sur¬ 
faces  in  steady  subsonic  flow  are  analysed.  The  methods  for  treating  the  singularities 
which  occur  in  the  kernel  of  linear  integral  equation  and  in  the  lift  distribution 
functions  are  also  discussed,  with  the  numerical  solution  of  the  integral  equation  given 
The  method  proposed  may  be  used  to  calculate  the  lift  distribution  of  the  wing  alone 
and  of  the  wing  with  full  or  partial  span  control  surfaces  on  both  the  leading  and 
the  trailing  edge.  The  numerical  results  are  in  good  agreement  with  experimental 
data,  and  are  as  accurate  as  those  obtuiaed  by  other  theories.  As  compared  with  the 
vortex-lattice  method,  both  the  required  computing  time  aud  the  computer  capacity 
are  reduced.  Thus,  the  numerical  calculations  may  be  carried  out  on  smaller  computers. 


50 


ON  THE  CAPACITY  OF  HOLOGRAPHIC  PHAGE-SHIFT  INTERFEROMETRIC 
TECHNIQUE  IN  THE  VISUALIZATION  OF  LOW  DENSITY  FLOWS* 


Li  Hua  yu,  Xu  Chao  yi ,  Shu  Jir.u,  Hu  Jin  min 
(Institute  of  Mechanics,  Academia  Sinica) 


ABSTRACT 


Gas  flows  at  high  Mach  numbers  are  usually  associate  with  very  low  gas  density, 
especially  in  simulating  flight  conditions  at  high  altitudes.  In  visualizing  such  flows 
the  question  usually  arises  as  to  how  the  sensitivity  of  the  optical  method  may  be 
increased. 

Paper  [2]  has  pointed  out  that,  if  the  reference  beam  of  holography  has  a  certain 
phase-shift  between  the  two  exposures,  a  higher  sensitivity  can  be  achieved.  Paper 
(2]  appends  an  experimental  result  in  which  the  phase-shift  value  is  jr/2.  The  pre¬ 
sent  paper  analyses  the  phusc-shift  interferometric  technique  in  detail.  According 
to  this  analysis,  when  the  phase-shift  value  is  tr/2,  a  resolution  limit  about  A/1000 
can  be  achieved.  This  is  25.4  times  higher  than  the  common  double-exposure  holog¬ 
raphic  method  in  which  the  phase-shift  value  is  zero.  Furthermore,  the  present  paper 
points  out  that  when  the  phase-shift  value  increases  from  x/2  to  n,  the  sensitivity 
also  increase  monotonously  (for  instance,  when  the  phase-shift  value  is  0.89  n, 'a  resolu¬ 
tion  limit  about  A/6000  can  be  achieved).  The  optimum  phase-shift  value  is  probably 
near  jt.  ...  • 

This  paper  presents  some  experimental  results  of  ft  low  density  flow.  Some  pro¬ 
blems  in  application  of  this  technique  ore  discussed. 


* 
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A  highly  sensitive  holographic  interferometric  method  deserv¬ 
ing  our  attention  is  introduced  in  this  paper.  As  pointed  out  in 
[2],  the  sensitivity  may  be  raised  when  there  is  a  phase  shift  in 
the  reference  beam  between  the  two  exposures  in  holographic  inter¬ 
ferometry.  The  experimental  results  for  a  phase  shift  of  v/2  is 
also  reported  in  that  paper.  In  this  paper,  we  analyze  in  detail 
the  interferometric  technique  and  point  out  that  at  a  chase  shift 
of  n/2,  this  method  may  achieve  a  resolution  of  about  1/1000  wave 
length,  25.^  times  the  sensitivity  of  ordinary  double  exposure 
holographic  technique  (with  0  phase  shift).  It  is  also  noted  that 
the  sensitivity  may  be  further  improved  for  phase  shift  greater 
than  tt/2  and  smaller  than  tt,  the  best  phase  shift  value  being  in 
the  neighborhood  of  m .  We  append  in  this  paper  the  experimental 
result  for  a  low  density  flow  to  demonstrate  the  ability  of  this 
method  in  visualizing  low  density  variable  flow  field  and  discuss 
some  of  the  problems  encountered  in  its  application. 

The  finite  fringe  interferometric  method  used  in  the  investi¬ 
gation  of  flow  fields  has  generally  a  resolution  of  about  1/20 
wave  length.  In  practice,  there  exists  a  great  number  of  situa¬ 
tions  where  the  state  of  matter  has  a  variation  much  less  than  1/20 
wave  length,  such  as  the  study  of  supersonic  air  flow,  all  kinds 
of  low  speed  flows,  small  cross-section  tube  flow,  etc.  Hence  the 
improvement  of  the  sensitivity  of  the  interferometric  technique 
has  always  been  an  important  research  topic.  The  phase  complement¬ 
ation  method  of  [1]  is  an  outstanding  example  of  works  in  this 
area.  In  that  paDer,  a  highly  sensitive  visualization  method  has 
been  developed  by  applying  the  principle  of  spatial  complex  filtered 
wave,  including  the  application  of  the  phase  shift  technique. 

In  contrast  to  the  principle  of  spatially  filtered  waves  with 
phase  shift  only  in  local  waves,  [2]  introduced  the  technique  of 
phase  shift  for  the  whole  wave  to  improve  sensitivity  by  changing 
the  phase  of  the  reference  beam  between  the  two  exposures  of  the 
hologram.  It  also  exhibited  the  pictures  taken  with  the  high 


sensitivity  obtained  after  shifting  the  phase  by  v/2.  [3]  applied 

similar  principles  in  non-holographi c  situations  to  realize  a  high 
sensitivity  in  ordinary  interferometric  methods  by  controlling  the 
initial  phase  difference  between  the  two  interfering  beams  to  be 
around  m/2. 

We  discussed  the  phase  shift  technique  in  the  holographic  inter¬ 
ferometric  method  in  an  attempt  to  explain  the  physical  interaction 
of  the  phase  shift  method  and  the  relationship  between  the  phase 
shift  value  and  the  sensitivity  of  the  measurement.  We  also  com¬ 
pare  the  methods  in  [2,3]  and  note  that  the  best  possible  phase 
shift  value  is  not  tt/2,  but  is  close  to  tt.  The  experimental 
results  are  included  in  the  appendix  with  a  discussion  of  the  prob¬ 
lems  that  exist  in  the  application. 

I.  ANALYSIS  OF  THE  PHASE  SHIFT  TECHNIQUE 

1 .  Effect  of  changing  the  phase  of  the  reference  beam  between 
the  two  exposures  during  flow  field  visualization  on  the 
double  exposure  hologram 

Let  the  object  wave  at  the  first  exposure  be  eiou.r>  ,  the 
reference  wave  be  ,'«■•»>  t  the  object  wave  at  the  second  exposure 
be  changed  to  J  the  reference  wave  be  where  <P 

is  the  phase  change  caused  by  the  flow  field  and  c  is  the  phase 
shift  value. 

Assume  that  the  recording  medium  records  the  intensity  linear¬ 
ly  [1],  then  according  to  the  general  theoretical  treatment,  the 
following  conclusion  may  be  arrived  at:  (1)  even  if  the  medium 
should  be  non-linear,  the  following  discussion  will  not  be  affected 
because  in  the  off-axis  holography,  the  non-linear  effect  only  pro¬ 
duces  higher  order  diffraction  images  with  large  diffraction  angles. 
It  is  separate  spatially  from  the  first  order  image  that  is  to  be 
used.  When  the  processed  hologram  is  illuminated  with  reference 

d3 


R(x,y) ,  the  two  objective  waves  emitted  are  respectively  ,««».»> 
and  <  c  =  0  is  the  situation  for  ordinary  double 

exposures.  This  conclusion  indicates  that  during  the  second  expo¬ 
sure,  the  increased  phase  c  of  the  reference  beam  is  equivalent 
to:  (1)  the  phase  c  subtracted  from  the  measured  objective  wave 

in  the  ordinary  double  exposure  holographic  interferometry;  (2) 
the  uniform  initial  phase  difference  c  in  the  whole  interfering 
plane  of  the  two  interfering  waves  in  non-holographic  interference. 

2 .  Intensity  distribution  of  low  density  flow  lnterferogram 
(flow  field  that  only  Induces  small  phase  changes) 


The  intensity  distribution  of  two  beam  lnterferogram  is 

/(*.>)*=  -j  (1  +  eoslip(x,  y)  —  f]}  (1) 

where  IQ  is  the  intensity  at  9>***  0  when  c  =  0,  i.e.}  the  background 
intensity  of  the  lnterferogram  when  c  =  0.  The  words  (x,y)  will  be 
omitted  in  the  rest  of  this  paper. 


Omitting  all  terms  above  sp*  under  the  condition  of  weak  phase 
changes,  we  get 

/(*,  y)  |(1  "*■  co*c)  +  ipiiac  —  cose] 

From  the  above  equation,  the  distribution  of  I(x,y)  when  <•  —  0,  */2, 
is  as  follows : 

K*. y)-'.(i--f-)  forf“°  ] 


O  +  <p)  for  c  "  nil 
f(*.y)  “  for 


t  m  P 


(3) 


t 


The  background  intensity  of  the  lnterferogram  may  be  expressed  as 

K*.  y)lf-o  ■■  (1  +  cose)  (.*0 

This  is  a  parameter  used  to  discuss  sensitivity  and  a  function  of  c 
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(5) 


/(»>  y)l*-o  “  /|  forf  —  O 

/(*»  tor  t  «■  w/2 

*  t 

/(*«  y)!*-o  “  0  for  f  * 

Besides,  we  define  the  intensity  enhancement  AI  by  the  following 
definition : 

A/  •-  /(*,  y)  ~  /(*,  y)|^  —  ^  (?  ““  f - co,f)  (6  ) 


3 .  Relationship  between  the  phase  shift  c  and  the  observed 
sensitivity 


The  flow  field  interferogram  contrast  is  defined  by 


<p  *in  c  - - jj-  co*  c 

1  +  cost 


This  function  is  monotonic  when  c  <  ir.  When  c  =  tt,  j^—oo  •  But 
due  to  the  effects  of  scattering  and  diffraction,  etc.,  there 
exists  randomly  light  so  that  will  never  be  zero  to  produce 

an  absolutely  dark  background.  Thus  practically  M  °°)  .  From 
the  contrast  angle,  naturally  the  sensitivity  is  higher  when  c  =  it  . 
But,  the  intensity  enhancement  should  also  be  a  practical  factor 
to  be  considered  because  when  9>  is  very  small,  AI  will  also  be 
very  small,  making  it  possibly  too  small  to  be  recorded.  In  the 
following  we  will  further  investigate  the  relationship  between 
AI  and  c . 


From  Equation  (6),  we  know  that  AI  has  a  very  large  value 
when  <•  —  — 2/<p)  (when  f  is  very  small,  this  condition  is  equi¬ 
valent  to  f  *»  w/2  ),  and  a  very  small  value  when  c  =  0  and  in  the 

vicinity  of  it.  The  difference  may  be  seen  from  the  following 


equation : 


A/  - 

4 

for  t  0 

A/--JV 

for  f  ■"  * 

A/  - 

2 

for  e  «•  n/2 

(8) 
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Summarizing  the  above  discussions,  we  know  that  if  there  is 
no  problem  of  recording  AI ,  the  best  value  of  the  phase  shift  is 
it.  But  if  “P  is  very  small,  the  recording  of  AI  will  be  limited 
by  the  facts  as  pointed  out  by  Equation  (8),  then  the  best  phase 
shift  practically  should  be  determined  by  a  compromise  between  the 
maximum  value  of  M  and  the  limiting  value  of  AI;  and  this  may  fall 
between  tt/2  and  it,  depending  on  the  actual  experimental  condi¬ 
tions.  With  the  improvement  of  technology,  the  best  phase  shift 
value  will  approach  it.  For  example,  as  pointed  out  in  [4],  the 
sensitivity  may  be  improved  by  another  1-2  orders  of  magnitude  if 
the  light  intensity  is  received  directly  with  photoelectric. 


4 .  Estimation  of  achievable  sensitivity  when  c  =  tt/2 


1)  Observable  sensitivity:  the  smallest  intensity  variation 
that  can  be  discerned  in  an  interferogram  by  directly  observing  it 
with  the  human  eye  is  called  the  observable  sensitivity.  It  Is 
generally  accepted  that  M  =  1/100  is  the  limit  discernible  by  the 
human  eye.  From  this  may  be  obtained  the  smallest  discernible 
phase  shift.  From  Equation  (7),  we  know 

«Ui-»  (9) 

Substituting  M  =  ±1/100  into  the  above,  we  get 


SP-'.U-  ^-f°r 
f  !<—■/>  ™  — ^for  | 


X 

31.4 

,  JL 

'  628 


where  Is  the  smallest  discernible  phase  change,  ALm^n  is  the 

corresponding  optical  path  differences,  X  Is  the  wave  length. 

From  this  it  can  be  seen  that  the  sensitivity  of  the  c  =  7r/2  method 
is  20  times  that  of  ordinary  double  exposure  method. 


2)  Recordable  sensitivity.  The  smallest  density  variation 
that  can  be  discerned  on  the  film  when  the  Interferogram  is  recorded 
photographically  is  called  the  recordable  sensitivity.  In  practice 


the  sensitivity  may  be  further  Improved  Lv  using  high  contrast 
film.  But  the  processing  requirement  is  more  demanding:  The 
exposure  must  be  correct,  the  control 
on  the  y  value  during  developing  must 
be  strict  to  guarantee  the  linearity 
of  the  record  and  the  necessary  value 
of  y •  Starting  with  the  characteris¬ 
tic  curve  of  the  recording  film  (Fig¬ 
ure  1),  we  have 


AD  •“  yllogC/Zo  +  AM)  — 

_o.,„rta(i  +  SH) 

logHd 

0.434r  — 

7  H» 

(10) 

Figure  1. 
tic  curve 
film 

CharacteriS' 
of  recording 

A H/H,  -  A/ 

(ID 

where  D  is  the  film  density,  H  is  the  film  exposure  value,  Hq  is 
the  film  background  exposure  value,  y  is  the  reverse  difference 
coefficient  of  the  film  after  developing,  AH  is  the  exposure  value 
enhancement  corresponding  to  AI,  and  AD  is  the  density  enhancement 
corresponding  to  AH.  In  Equation  (10),  the  condition  for  the 
approximation  is  that  AH/Hq  is  very  small.  Solving  Equations  (9)- 
(11)  simultaneously,  we  get 

ADU - o.43ir-^-.  ADI-*,- 0.431  r*  (12) 

Under  the  condition  that  AD  |  <3iscemible  =  1/100  ^  -  and  Y  =  11 »  we  get 


9nu|r«« 


2n 

41.3’ 


X 

41.3 


Thus,  the  sensitivity  of  the  c  =  m/2  method  is  26.4  times  that 
of  the  ordinary  double  exposure  method.  In  this  discussion,  the 
reason  why  the  recordable  sensitivity  Is  higher  than  the  observable 
sensitivity  is  that  films  with  high  y  value  have  been  used  (contrast 
Equations  (9)  and  (12)).  If  r  <  1/0.434  —  2.31  ,  then  the  record¬ 

able  sensitivity  will  be  lower  than  the  observable  sensitivity. 


D=log  0,  0=F  /F  is  called  the  optical  resistivity.  Fq 
dent  light  transmittance  on  the  film  when  measuring  D. 
optical  transmittance  through  the  film.  AD=log(0/0  ). 
1/1000,  0/0o-1.023. 


is  the  inci- 
F  is  the 
When  AD= 
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y"-1*  for  ordinary  films,  therefore,  ordinarily  the  recordable 
sensitivity  is  considered  to  be  lower. 

The  above  method  was  also  used  in  [3]  to  estimate  the  record¬ 
able  sensitivity  of  its  optical  system.  According  to  its  estimate, 
the  highest  discernible  value  for  its  system  is  A/500,  exactly  half 
of  the  value  in  our  paper.  The  reason  is  that  complementary  double 
exposure  was  used  in  the  system  to  compensate  for  the  original 
error  in  the  optical  elements  so  that  the  background  intensity 
could  be  made  uniform.  Thus,  because  the  phase  shift  was  fixed 
near  n/2,  AH  was  similar  to  that  in  our  paper  but  Hq  was  doubled. 
Hence  the  highest  discernible  value  is  lowered  by  one  half.  Fur¬ 
thermore,  holographic  interferometric  technique  can  maintain  the 
same  phase  shift  for  every  point  in  the  flow  field.  The  method  in 
[3]  could  not  achieve  this  due  to  error  in  the  processing.  There¬ 
fore,  the  sensitivity  at  each  point  of  the  flow  field  is  different. 
The  holographic  method  is  also  superior  in  this  regard. 

In  theory,  the  sensitivity  at  phase  shift  tt/2  is  lower  than 
that  at  phase  shift  ir.  But  there  are  two  advantages  for  phase 
shift  tt/2  :  1)  AI  is  largest,  making  it  easier  to  record;  2)  AI 

and  <P  are  linearly  related  eo  that  quantitative  analysis  is 
easier.  Since  its  sensitivity  is  already  high,  this  will  Improve 
its  practical  value. 

5 .  Estimation  of  the  sensitivity  when  c  approaches  it 

Since  we  have  only  considered  the  effects  of  the  major  fac¬ 
tors  in  our  theory,  we  have  reached  the  conclusion  that  M«»oo 
when  c  =  tt.  In  reality,  the  closer  is  c  to  it,  the  greater  are  the 
effects  of  those  factors  not  included  in  the  theory.  For  example, 
when  W*>°o  t  infinitesimally  small  AI  can  also  be  recorded.  Here 
the  exposure  time  should  be  infinitely  long,  but  the  fogginess  of 
the  film,  the  stray  light  during  reconstruction  and  other  factors 
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will  make  the  recording  meaningless.  That  is  to  say,  if  the  V> 
that  needs  to  be  recorded  is  very  small,  then  the  best  phase  shift 
value  must  be  smaller  than  n  and  must  be  a  compromise  between  the 
limiting  value  of  AI  and  the  largest  value  of  M.  Two  examples 
will  be  used  below  to  estimate  the  recordable  sensitivity. 


1)  The  recordable  sensitivity  when  c  =  150°. 
ing  y  «•  AD  —  l/ioo  into  Equations  (7),  (10)  and  (11) 


2% 

5036' 


5036 


By  substitut- 
,  we  get 


The  intensity  enhancement  is  very  small  for  such  a  small  phase 
shift.  We  compare  them  as  follows: 


A/ 1  f  Ml. «.*•>!«/»*  “  7.784  X  10-’  X 
AI 1 “  1.248  X  1 0"’  X 

2 

Al  I  r—IXT-r—ln/Xm  “  6.254  X  10-*  X 


2)  Recordable  sensitivity  when  c  =  160°.  Similarly,  the 
following  data  are  obtained 


fall  I  «m|«0* 


2* 

6195* 


AI  l,M>...r-l«/*lM  “ 

AI  | 

AI  | 


At.,.U-.-6195 

5.14  x  10-7  x 

2 

1.014  X  10"7  X 


3  473  X  10-4  X 


It  can,  therefore,  be  seen  that  the  sensitivity  continues  to 
improve  and  AI  to  decrease  when  the  phase  shift  approaches  tt  but 
when  the  phase  shift  deviates  slightly  from  it,  the  value  of  AI 
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will  be  increased  greatly  over  its  value  at  c  =  u,  thus  reaching 
a  stage  where  it  may  be  usable. 

II.  RESULTS  AND  EXPERIMENTAL  SET-UP 

1 .  Experiment  set-up 


The  following  experiment  was  performed  to  investigate  the 
feasibility  of  applying  the  phase  shift  interferometric  technique 
in  the  visualization  of  low  density  flow  field.  The  object  of  the 
experiment  is  a  small  column  of  freely  flowing  air,  passing  through 
a  tube  under  constant  total  pressure  of  700  mm  Hg  in  a  container 
into  a  cylindrical  exhaust  tube  of  diameter  1  mm  and  length  4  mm 
and  then  into  the  atmosphere.  The  commonest  off-axis  holygraphic 
optical  path  Is  used  (Figure  2).  The  phase  shift  is  achieved  by 
using  a  phase  regulator  in  the  path  of  the  reference  beam  prior  to 
b earn- expans ion 


Figure  2.  Diagram  of  experimental  optical  path. 

(1)  HeNe  laser  beam;  H  holographic  plate;  S  beam  divider;  M1M2 
total  reflecting  mirrors;  LiL2  beam  expanding  lenses;  L3  parallel 
focusing  lens;  P  phase  regulator;  E  exhaust  tube 


The  phase  regulator  is  a  </>12mm  x  100mm  metal  cavity  with 
the  two  ends  sealed  by  optical  glass.  The  pressure  in  the  cavity 
is  very  stably  controlled  by  a  two  stage  pressure  stabilizer. 


The  various  phase  shift  values  c  =  0,  ±tt/2,  *  are  achieved  con¬ 
veniently  by  controlling  the  pressure  value.  The  sign  of  the 
phase  shift  Is  determined  by  filling  the  phase  regulator  with  air 
either  during  the  first  or  the  second  exposure.  The  distortion 
of  the  optical  glasses  at  the  two  ends  of  the  phase  regulator 
should  be  minimized  in  the  design  of  the  regulator.  This  problem 
is  solved  in  our  case  by  controlling  the  length  of  the  metal 
cavity.  Two  factors  are  taken  into  consideration:  1)  the  applied 
pressure  should  be  minimized.  A  pressure  of  only  60  mm  HgO  is 
required  to  achieve  a  phase  shift  of  u/2  in  our  regulator  (at  room 
temperature  20°C)  1^;  2)  the  air  pressure  should  be  within  the 

operating  range  of  the  pressure  stabilizer.  Therefore,  the  cavity 
cannot  be  too  long.  There  are  other  schemes  for  the  phase  regulator. 
For  example,  the  piezoelectric  crystal  is  also  a  possibility. 


2.  Experimental  results 


Experiment  indicates  that  the  observed  sensitivity  can  really 
be  greatly  improved  with  the  phase  shift  technique.  Limited  by 
circumstances,  we  recorded  with  ordinary  aeronautic  microfilm  with 
no  control  on  the  y  value.  The  recordable  sensitivity  is  lower 
than  the  observable  sensitivity  and  far  lower  than  the  standard 
in  the  above  discussion.  For  results,  please  refer  to  the  pictures 
1-6  in  Plate  I.  Picture  1:  Recorded  with  the  ordinary  double 
exposure  method.  A  very  short  dark  flow  column  could  be  vaguely 
discerned  during  observation.  Nothing  can  be  seen  in  the  photo¬ 
graphic  record.  Picture  2:  Recorded  by  the  finite  fringe  double 
exposure  method.  The  limited  fringes  are  produced  by  the  liquid 
wedge.  A  disturbance  with  an  optical  path  difference  of  about 
A/10  can  be  seen  at  the  exit  of  the  flow  column.  The  bending  at 


^The  pressure  required  for  phase  shift  tt/2  may  be  calculated  with 
the  following  formula:  t  lft  T 

P  Tf  Vi  “ 


where  K=2.9*4  x  10~\  L  is  the  length  of  the  cavity,  T  is  abso¬ 
lute  temperature,  p  ,  T  are  the  pressure  and  absolute  tempera¬ 
ture  at  STP .  After°subst ituting  in  the  values,  we  have 
o  *  0.0201<T,  dimension  being  mm  H«0. 

v  relative  2 
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the  second  fringe  is  already  very  small.  Pictures  3  and  4 : 

Recorded  by  the  double  exposure  method  with  phase  shift  ±tt/2, 

A  white  (corresponding  to  a  phase  shift  ti/2)  and  a  dark  (corres¬ 
ponding  to  a  phase  shift  of  n/2)  flow  column  can  be  individually 
observed  clearly  within  a  fairly  long  range.  The  visualization 
capability  is  much  higher  than  that  in  pictures  1  and  2.  Picture 
5:  Recorded  by  the  double  exposure  method  with  phase  shift  r. 

A  white  flow  column  over  a  dark  blackground  can  be  seen.  Picture 
6:  The  experiment  was  disturbed  by  some  unusual  factors,  but 

they  accidentally  make  the  phase  shift  in  the  flow  region  come 
out  to  be  ideally  tt.  The  resulting  visualization  of  the  flow  is 
the  best. 

We  also  tested  the  operation  of  the  phase  regulator  in  real 
time.  When  the  pressure  is  increased  continuously,  the  fringes 
continue  to  scan  with  the  value  of  position  shift  in  agreement 
with  calculated  value.  Therefore,  it  has  operated  normally. 

3 .  Discussion  of  existing  problems 

In  practical  applications,  a  serious  impediment  of  our  method 
is  the  stringent  requirement  on  the  environment .  The  experiment  in 
our  paper  had  been  performed  over  80  times  with  satisfactory  results 
In  only  a  few  cases.  The  reason  is  that  the  disturbance  of  the 
surrounding  can  easily  exceed  tt/2  (e.g.,  air  flow,  thermal  dis¬ 
turbance,  small  vibrations,...)  so  that  the  desired  phase  shift 
value  cannot  be  achieved,  leading  to  the  failure  of  the  experi¬ 
ment  or  to  disturbances  of  various  degree.  If  we  use  the  single 
exposure  hologram  as  a  wave  memory  storage  and  then  use  it  for 
image  reconstruction  with  the  phase  shift  method  in  a  better  en¬ 
vironment  may  be  a  successful  way  to  do  this  experiment.  This 
technique  has  a  good  prospect.  Its  development  awaits  future 
practices . 
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The  second  problem  is  that  there  still  exists  certain  diffi¬ 
culties  in  using  this  technique  in  quantitative  measurement.  But 
if  the  error  can  be  reduced,  then  by  usinn  the  linear  relation¬ 
ship  between  intensity  and  phase,  the  quantitative  measurement 
may  be  achieved  with  a  micro-densitometer. 

IV.  CONCLUSIONS 

The  special  feature  of  the  method  in  this  paper  is  its  high 
sensitivity  and  theoretical  simplicity.  Concerning  the  latter, 
we  have  not  seen  any  paper  discussing  it,  probably  because  that  it 
is  difficult  to  achieve  in  practice.  Theoretically,  the  difficulty 
to  achieve  constant  phase  shift  has  been  eliminated  with  the  appear¬ 
ance  of  holographic  photography.  Hence,  it  can  be  predicted  that 
the  high  sensitivity  of  our  method  will  find  Its  use  in  a  suitable 
field. 
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ANALYTICAL  DESIGN  FOR  INTERNAL  BURNING  STAR  GRAINS  OF  SOLID  ROCKETS* 


Lu  Changtang 

For  the  design  of  the  internal  burning  star  grains  of  solid 
rocket  motors,  in  the  past  people  usually  carried  out  pure  geo¬ 
metric  research  and  the  corresponding  trial  and  error  method. 

The  formulae  and  computational  curves  in  [1]  are  typical.  They 
are  widely  used  in  the  engineering  design  and  scientific  educa¬ 
tion  in  foreign  countries.  Practical  experience  shows  that  the 
pure  geometric  formulae  and  curves  in  [1]  are  complex  and  numerous 
and  yet  incomplete;  the  corresponding  trial  and  error  method  not 
only  involves  large  volumes  of  computation,  but  is  also  uncertain, 
difficult  to  guarantee  the  accuracy  of  various  specifications  and 
to  achieve,  in  particular,  the  optimal  design  standard.  In  this 
paper,  we  try  to  use  analytic  design  instead  of  trial  and  error. 

Different  from  the  traditional  pure  geo¬ 
metric  research,  we  organically  combine 
the  various  geometric  parameters  of  the 
Internal  burning  star  grain,  the  various 
characteristic  parameters  of  solid  pro¬ 
pellant,  and  the  technical  specifications 
of  the  rocket  motor  together,  establish  the 
set  of  design  equations  according  to  the 
best  principle  of  grain  design,  and  then 
solve  the  set  of  equations  quickly  by  in¬ 
corporating  simple  equations  and  curves. 
Compared  with  [1-3] ,  our  method  not  only 
involves  less  computational  work  and  guar¬ 
antees  the  required  technical  specifica¬ 
tions  but  also  yields  a  design  that 
approaches  the  optimal,  improves  the  char¬ 
acteristics  of  the  rocket  motor,  saves 


Figure  1.  Diagram  of 
design  unit  for  the 
internal  burning  star 
grain 


I - 

Received  April  25,  1978 


64 


materials  and  requires  less  difficult  work  on  the  computer. 

1 .  The  ,-*oor.'" 1  rlc  functional  relations  of  the  Internal  burning 
star  grain 


As  in  Figure  1,  it  is  not  difficult  to  obtain  the  functional 
relationship  between  the  geometric  parameters  of  the  internal 
burning  star  grain  according  to  Piobert's  combustion  law  as  done 
in  [1].  For  simplicity  in  the  calculation,  we  first  give  below 
only  the  function  k  (Figure  2)  and  function  sq  (for  the  meaning  of 
symbols,  see  Figure  1)  which  are  only  related  to  N,  the  number  of 
points  in  the  star  and  9/2,  the  initial  half  angle  of  the  star 
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From  these,  the  combustion  surface  in  the  first  stage  and  the 

initial  port  area  Apol  may  be  simplified  into  the  two  equations 
below : 

Atl  -  L,U  +  */.  +  2 *(y  -  /,))  C3) 

An”sJ  +  Pc  +  nfi  -  0.5*(/'  +  ID  (4 ) 


where  is  the  grain  length.  For  other  symbols,  see  Figure  1, 
For  value  of  c,  see  the  equation  below: 

f  ”  ( i  —  e  )ir  d-  N  sin  e  —  cos  e  — ■  —  N  sin's  —  ctg  —  (-5  ) 

N  N  N  2 


From  analysis,  we  know  that  angle  8/2,  called  the  dynamic 
angular  variable,  increases  monotonically  in  the  interval  [0J2, 0„./2] 
with  the  combustion  time,  starting  from  y@o/2,  the  star  edge  van¬ 
ishing  point.  0/2  is  the  star  half  angle  at  the  end  of  combustion. 

W 

0/2  is  related  to  the  time-varying  y  as  follows: 


—  —  arecos 


(y  +  /)/< 


(6) 
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By  usinm  the  dynamic  angular  variable  0/2,  th^  r<-^^nd  stare 
combustion  area  A^-q  and  the  remaining  grain  area  may  be 
simplified  into  the  following  functional  forms: 


-*•”  - 1.  [*-  +  w  (f  -  f )]  '-^9>  +  ML'(‘  - 

Ai  —  > «/(D,  —  /)e  +  iVClt'  4-  /)'  (  —  —  —  —  *in  -^-cos  — 1 —  NF  sin  t  —  co*  t  — 
'  \  2  2  2  2/  y  N 


(7) 

(8) 


where  the  remaining  grain  surface  star  half  angle  0w/2  at  the  end 
of  combustion  is 
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It  can  be  shown  from  mathematical  analysis  that  Abll  has  a 

minimum  value  A.  at  Q/2  and  is  infinite  stiffness  form  (Figure 
bmin  “ 

3  is  the  functional  curve  of  its  infinite  stiffness  form)  is 

n 

llo  £  77 

Jsia--2N - rr~r  +  2*(1  —  e)  (1°) 

/  .in  (tf/2) 

Similarly,  it  may  be  proved  that  A^  has  a  minimum  value  Afmin 

at  Ox.~JL-.JL  i.e.  ,  W  +  j  ^  . 

2  2  N  I 

The  existence  of  A^min  and  Afmin  is  an  important  feature  of 
internal  burning  star  grain. 


The  effective  grain  area  Acq  when  there  is  no  star  circular 
angle,  namely  f ^  =  0 ,  is  very  important  to  the  establishment  and 
solution  of  the  design  equations.  The  following,  formula  gives  its 
value : 

->  t  oi "  +  0  ~  °  'n  ~1'w"  ~  ””‘5 

+  Nl1  ctg  un'c  -  N(.W  +  /)>  •  —  y  —  *in  y1  co* 
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\jbt m**  -  » **“ 


when  there  is  star  circular 


Then  the  effective  grain  area  A 
ancle,  namely  f^  \  0  is 

AM  -  A*  -  0.5(2«  -  0/i  (13) 


The  minimum  radius  Rq  of  the  internal  burning  star  grain 
when  f^  =  0  gives  the  governing  relationship  between  the  geome¬ 
tric  parameters  of  the  star  grain.  Since  it  is  always  true  that 
^Jidn<0  ,  therefore,  R  is  a  decreasing  function  of  the 

angular  coefficient  e.  Let  Rq  >  0 ,  then  the  upper 

limit  eM  of  £  is: 

Nap 


?(■ 


(1*0 


Whenever  e<eNt  ,  then  R  >0  and  R  ,  >  0  (when  fn  k  0). 

o  ol  1 


2.  Establishment  of  the  set  of  main  design  equations 


The  known  conditions  are: 

1)  The  required  technical  specifications  for  the  motor: 

The  total  impuse  or  required  range,  law  of  propulsion  Fmin_ 

Fmax  operating  time  tmln-tmax,  operating  temperature  Tmin-Tmax°C 
and  the  outer  diameter  D  of  the  rocket  motor  or  its  limit,  etc. 

2)  Property  of  solid  propellant:  specific  impulse  I  ,  com- 

n 

bustion  rate  y  =  ap  ,  temperature  coefficient  a  ,  limit  of  air 
current  sensitivity  coefficient  k^ ,  propulsion  coefficient  C^,, 
characteristic  speed  c  ,  critical  pressure  ,  specific  gravity 
p,  etc. 

The  optimal  grain  design  principle  obeyed  by  the  set  of 
design  equations:  to  be  able  to  guarantee  the  relization  of  the 
total  impulse,  propulsion  power  requirement,  operating  work,  etc; 
to  be  able  to  burn  steadily  and  normally;  to  have  a  large  effective 
grain  coefficient  n  and  a  small  coefficient  nf,  etc.  The  stresses 

v  1 
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under  conditions  permitted  by  the  grain  stress. 

Thus,  the  design  equation  that  guarantees  the  requirement 
total  impulse  of  the  rocket  motor  is 

l*f[Alt  —  0.5(2*  —  ^)/f]  ™  V,  (17) 

The  design  equation  guaranteeing  the  value  of  the  lower  limit  of 
propulsion  power  is 


2  N 


tin  e  t: 

— f  +  C1-8^ 


*4 ♦ml* 

\c; 


(18) 


The  design  equation  guaranteeing  the  initial  peak  pressure 
limit  for  safe,  stable  combustion  is 


U-Qg-QML.  _ 
A „  +  0.5(2»  -  O/l 


(19) 


By  eliminating  the  grain  length  Lp  and  f1  and  substituting 
into  the  above  Equations  (1),  (2),  (.8)  and  (12),  the  equation  of 
the  angular  coefficient  e  may  be  obtained 


fcpD'  -  ^.,,/eCD,  -  /)  -  Mt.pOk  +  \y 


arccos- 


sr 

I  ,'n877 


W  +f  2 


/*inev  /  /sinerjA  _  • 

— - —  (inlarccot - J  +  N<c«pB  sine  — cose - e*.nV,Ai£i.2xKl  —  e) 

W  +  /  \  W  +  f  /J  ^  N  N 

w 

r  dnB  at 

•  [w/^— f 


-2  NA&JV.k.p- 


.  d 

sin  — 
2 


n 

sins  T7 

2M  - -  +  2xf(l  -  e ) 


.  fl 

tin  — 
2 


C2Q) 


+  2*M,..u(l-0  +  2NM».U  +  4^1 

\2  2  “  2/J  * 

JD\  -  WKV  +  /) 

V. 

-  2Nctg  -  8jr,/lks  -  ?-fl-  sins—  4-  4*V/*ctg  ^ 
2/  Jail  W  2 


+  A*, 


-  2  nf>  ( 


sr 

sine  tt 

2N/  - £-  +  2*/(l  -  e) 

sin- 

2 


Ahr  +  2*  -  2N 
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X  tin’s 


X  /  M\  * 

_  /tinSTT  I  /tinSTTl/tinSTT 

- +  /)’  arcco* - —  —  —  tlnltrccos - 1 - 

N  L  W  +  j  2  \  W  ]  W  +  t 


n 

tins  tt 


-  ta  2M  +  2*/(l  -  «)  -y ND>'  *  (f  +  Jj  -  f  ~  «*  f ) 

•in  — 

2  4  . 

+  1.W1I  +  I)  ■  (i  +  i  -  |  -  «l  f)  +  IN'/*  (f  +  j  -  f  -  «*  f)’ 

+  \hNIW 6  -  Ctg  —)  +  4/V’//  (—  +  —  -  -  etj-S) 

V  2  N  2  2/  \2  N  Z  *  2  / 

* 

x  _  2/v’/»  (i  +  ~  -  &  -  ctg  ^  ctg  .  tin’s  —  +  2N‘ (W  +  /)’ 

.  0,  \2  N  2  8  2  /  8  2  -  N 


li,neJj  ,  .  (  /«ne^V-^/x  ,  *  0,  0,\ 

- tin  \  arc  cos - ) -  1  —  + - *  —  ctg  — *  ) 

I  V+f  2  \  W+l/w+j\\Z  N  2  2> 


X  arccos 


+  2NV'A&>  *  (“  +  T,  ~  4"  “  ctfi  T4)  2*/(l  -  e).+ 

'2  N  2  2  / 


tine  tv 

2 Nl - -  •  -  0(20) 

.  0 
tin  — 

2 


It  may  be  seen  that  this  equation  contains  simultaneously  all 
the  factors  such  as  the  required  technical  specifications,  the 
property  of  the  propellant  and  genetric  parameters  of  the  inter¬ 
nal  burning  grain,  etc. 

3 .  Requirement  of  the  design  equations  and  obtaining  the 
geometric  parameter  of  the  grain 

Obviously,  it  is  difficult  to  solve  for  e  In  Equation  (20). 

We  change  it  into  two  functional  curves  Y-^,  Y£  with  e  as  the 

variable:  _ 

Y, O’, Yj-L,[s,-v/2A(2x-01  (21) 


where 


The  necessary  equation  with  A  as  solution  can  only  be  solved  when 
A  >  0. 


For  each  positive  integer  N,  the  angle  9Q/2  has  an  upper  limit 
0»ri/2  that  makes  A  =  0  and  f,  —  O,0wk/2  ,  is  an  increasing  function 
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of  N,  therefore,  the  angle  0/2  must  be  within  the  interval 
(0,  and  be  determined  according  to  the  initial  power 

requirement.  N  and  0q/2  are  called  the  parametric  variables  of 
Equation  (20).  To  obtain  good  overall  grain  property,  the 
values  of  N  and  9q/2  should  be  used  (such  as  A/ —  2,l,dj2  =  10°). 

After  N  and  6q/2  are  fixed,  use  (N,  0q/2)  to  search  k  value  from 
Figure  2  and  then  calculate  the  curves  Y  ,  Y^  of  e:  For  each  e 
selected,  find  the  value  of  smin/l  from  Figure  3  with  (N,  e)  and 
substitute  into  the  following  equation  to  get  Lp : 

L,  -  -dlKla- 

l.isla'  (23) 

-  / 


Figure  2.  Diagram  of  functions  to  compute  k. 

From  this  we  may  obtain  V  /L  and  then  substitute  e  and  N, 

c  p 

9/2,  k  into  Equation  12  to  get  A__  and  hence  A.  When  A  >  0,  the 

O  CO 

functions  Y1  and  Y^  corresponding  to  e  may  be  obtained  from  Equa¬ 
tion  (21)  and  s„  and  A,  which  are  obtained  for  Equations  (2)  and 
(8).  A  series  of  functions  Y-^  and  Yg  of  e  may  be  obtained  in  this 
way.  The  solution  e  solution  of  Equation  (20)  is  given  by  the 
abscissa  of  the  point  of  intersection  of  the  decreasing  function 
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Y1  and  increasing  function  Y2  of  e.  The  angular  coefficient  of 
the  grain  can  then  be  obtained;  the  vertical  coordinate  is  the 


Ab01  of  the  r'rain-  sr.in//1  is  found  from  FJCure  3  by  using  (N, 

£  ,  ^  ).  As  before,  we  may  calculate  the  grain  length  L  and 

solution  P 


,  etc.  Vc/L  is  the  value  of  Ac01. 


Figure  3.  Diagram  of  the  function  s  .  /I, 


We  then  calculate  f  and  R  .  from  the  following  two  equations: 

1  ot 


R«  —  c*c-* 


When  f^  =  0,  we  know  from  Equation  (.1*0  that  whenever  £<£Nyp 

then  certainly  Rq  >  0,  and  when  f ^  K  0 ,  then  if  the  equation  has 

a  solution,  R  ^  >  0.  Furthermore,  the  smaller  is  90/2,  the  larger 

is  R  .. 
ol 

Discussions  on  the  solution: 


When  we  vary  e  and  A  <  0  is  still  true,  then  the  equation  has 
no  solution  and  we  must  change  to  larger  N  and  smaller  90/2 .  We 


can  also  vary  the  known  parameters  such  as  k^  within  allowed 

limits  (Notice:  increases  when  N  is  large.  A^  decreases 

when  0Q/2  is  large.  When  necessary,  v/e  may  adjust  values  of  N  and 

0/2 ) . 
o 

By  now  all  the  geometric  parameters  of  internal  burning  grain 
has  been  found.  Finally,  for  the  law  governing  the  variation  of 
the  combustion  surface: 


(0,  y»yi]  is  the  first  stage  of  combustion,  i.e.,  the  linear 


stage 


■-  •  r’  *•  * 

ltla,N 

co#~» 

2 


In  the  interval  [0,^],  we  have 

Ah  li*  +  Vi  +  2«<y  —  /,)]!,  (26) 

In  the  interval  (A,  y*jt } .  we  have 

A„-0>  +  Ky)L,  (27) 

lyy>.  M'l  is  the  second  stage  of  combustion,  i.e.,  the  non-linear 
stage.  The  variation  of  A^-^  may  be  found  readily  from  Equation 
(7)  using  the  dynamic  angular  variable  6/2. 


It  should  be  pointed  out,  in  passing,  that  the  statement  by 
F.  A.  Williams  et  al  in  [2]:  "In  the  second  stage  of  combustion, 
...the  grain  is  ,  i.e.,  A^  increases  linearly  according 

to  y  =  rt"  is  not  true.  During  the  second  stage  A^-q  is  non¬ 
linear  and  when  6^/2  <  6/2  it  is  monotonically  decreasing  at  first 
and  then  monotonically  increasing. 
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NUMERICAL  COMPUTATION  AND  ANALYSIS  OF  THE  FLOW  FIELD  IN  A  LARGE 
SHOCK  TUBE  WITH  A  VARIABLE  CROSS-SECTIONAL  AREA* 


Li  Wen  xuan,  Wane  Jiajun 
(Beijing  University) 

Reference  [l]  has  clarified  and  explained  certain  flow  char¬ 
acteristics  of  the  new  shock  wave  that  exists  in  shock  tubes  with 
violently  varying  cross-section.  These  characteristics  are  what 
the  approximate  analysis  of  Chisnell,  Chester,  Whitman  and  others 
cannot  describe  because  they  have  restricted  that  the  cross-sec¬ 
tion  changes  slowly  and  that  there  should  not  be  many  strong  dis¬ 
continuities.  This  paper  has  taken  into  consideration  the  state 
equation  of  the  explosive  gas  produced  by  the  combustion  of  the 
explosive  in  the  high  pressure  section  as  well  as  the  contact  dis¬ 
continuity  between  the  explosive  gas  and  air  and  also  extended  the 
calculation  to  the  whole  tube  (Figure  1).  In  order  to  treat  many 
strong  discontinuities,  we  adopted  the  difference  method  with  arti¬ 
ficial  viscosity  term  q.  But  different  from  [l],  the  Lagrangian 
coordinates  are  used  In  this  paper  so  that  the  contact  disconti¬ 
nuities  can  be  easily  treated.  Also  the  head  of  the  major  shock 
wave  has  little  vibrations,  occupying  only  2-3  grid  spacings. 


T 

Ti 

• 

U.,  -  J 

*  tttSBMk  •- 
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Figure  1.  Diagram  of  the  shock  tube 
Key:  1 — high  pressure  section;  2 — film;  3 — medium  pressure  section; 
*•— •ection  with  variable  cross-sectional  area;  5 — transitional 
section;  6 — experimental  section;  7 — exit  section; 8 — explosive  gas  in 
#high  pressure  section ; 9--gas  in  medium  and  low  pressure  sections 
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1.  DIKENGI0NLES3  EQUATIONS  AND  FINITE  DIFFERENCE  GRID 


The  dir.ensionless  equations  are 

dr  -  u-  &L.  *  £sAx.  in.  _ _ *  dP-  ,,  _  . 

di  ’  dR  PA  ’  di  .  .  *4  dK'.ft  ,.a.^>d<  (1.1; 

where  p  ,  pQ  are  respectively  the  static  air  pressure  and  the 
density;  d  is  the  diameter  of  the  thin  tube;  r  is  the  Euler  coor¬ 
dinate;  R  is  the  Lagrangian  coordinate;  t  is  the  time;  p,p,u  are 
respectively  the  pressure,  density  and  the  velocity  along  the  r 
direction;  e  is  the  internal  energy;  A(r)  is  the  tube  cross-sec¬ 
tional  area.  To  make  the  quantities  dimensionless,  they  have 
been  divided  respectively  by  d,  d,  d/y/JJ^,,  p,,  VpJp>.  h/p*.  d\ 

The  subscript  0  denotes  the  initial  values  of  the  various  quanti¬ 
ties.  The  difference  equations  are 


r?«‘  -r1_  .jy+l.  'In  ~.dll  ,  *±id lit 

I  "I  *  n  _  D 


At  Rj+i  ~~  Rt  • 


A t  Ri+i  Ri- i 


-»+l  _  . 

w.tvj 


c*i  -  K<i-  °::p 

...  r  a  < 0  w 

m‘-|0,  h  *<«#>.» 


where 

(«->;;}  - -  -r ?  *  c*>r  - 
-  C}-C} 
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2  2 
c  is  a  constant  between  1-4.  We  take  c 


2. 


Treatment  of  the  initial  values  "Y* 


Substituting  Equation 


u'r  - 


A  t_  -dL  (ikLUlzL 


(1.1)  as 
+  0(Al*); 


,  we  get 

di 

the  accuracy  of  „y>  is 


0(A/J). 


The  surface  separating  the  two  gases  is  treated  with  the 
method  of  [3].  The  advantage  of  using  the  Lagrangian  coordinates 
in  the  calculation  is  that  the  surface  separating  the  two  gases 
always  lies  at  the  grid  point  J,  thus  simplifying  the  treatment 
of  the  surface.  At  this  point,  there  is  discontinuity  in  p,e. 
Hence,  we  introduce  j_  /dp\  _  3(p)+t  -  tf-f)  -  y  -  Pj-0 

where  AR-^AR^  are  p,\dRJj  p.Atf,  +  p<AR,  * 

respectively  the  step  lengths  on  the  opposite  sides  of  the  separat¬ 
ing  surface;  are  respectively  the  initial  densities  of  the 

explosive  gas  and  air. 


Two  schemes  have  been  used  for  the  state  equation  of  the 
explosive  gas.  In  the  first  state  equation,  the  collisions  between 
hard  sphere  molecular  model  of  the  gas  is  taken  into  consideration. 
Then  we  take  a  state  equation  [2]  with  two  Wylie  coefficients.  We 
also  adopt  from  [4]  the  variable  thermodynamic  functions:  speci¬ 
fic  heat  at  constant  volume  <•,( T)  ,  specific  heat  at  constant 
pressure  *,(T)  >  speed  of  sound  tl(pi 7)  ,  internal  energy  e(T), 

etc.,  but  we  neglect  the  chemical  reactions  that  the  explosive  gas 
further  undergoes  after  t  >  0.  The  approximation  formulae  (T  being 
the  temperature)  for  each  section  of  the  equation  of  state  are 


*(7) 


p  -  1.108p7(l  +  2.5538  X  10'V  +  0.4045  X  KTV); 

.4.0808(7  -  3.4344)  +  0.1361(7  -  3.4344)*,  3.4344  <  7  '*  «  \ 

7.4089  +  4.5314(7—5.1516)  +  0.0846(7 -3.1516)’.  5.1516  <  7  <6.8688 

15.4399  +  4.8219(7-6.8688)  +  0.0543(7  -  6.8688)*,  6.8688  <  7  <  8.5860 
23.8799  +  5.0086(7-8.5860)  +  0.0377(7-8.5860)*,  8.5860  <  7  <  10.303 
'32^9517  +  5.1381(7 — 10.303)  +  0.02861(7-10.303)*,  10.303  <  7 
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<-.cn 


,3.6822  +  0.2456(  T  -  3.4344), 
4.0890  +  0.1524(7  -  5.1516), 
4.3508  +  0.09H0(T  —  6.8688), 
4.5193  +  0.0681(T  -  8.5860), 
4.6372  +  0.0516(T  —  10.303), 


T  <5.1516 
5.1516  <  T  <  6.8688 
6.8688  <  T  <  8.5860 
8.5860  <  T  <  10.303 
10.303  <  T 


a‘(p,  T)  —  1.108T 


[l  +  -=-(»  - 

l  e. 


0.2476  X  10 


V)] 


X  (1  +  5.1077  X  10- V>  +  1.2137  X  Mr*p») 


The  second  kind  of  state  equation  uses  the  perfect  gas  model 
with  equivalent  constant  specific  heat.  We  take  R  =  3.193J/kg, 
and  the  absolute  thermodynamic  index  y  =  ~  1.225  .  After  made 

dimensionless,  they  become  p  —  l.l08pT;  e  —  4.925T;  R  —  1.108. 

More  factors  are  considered  in  the  scheme  of  the  state  equa¬ 
tion  of  the  first  kind,  but  more  calculations  are  involved;  less 
factors  are  taken  into  consideration  for  the  second  scheme,  but 
the  calculations  involved  are  simple.  The  results  from  the  two 
schemes  are  close.  Especially  during  the  later  stages  of  the  motion 
the  results  basically  overlap.  In  the  complete  computational  pro¬ 
cess,  the  largest  difference  between  the  results  of  the  two  schemes 
only  appears  In  the  Initial  stage,  and  never  exceeds  5%  (Figure  4). 
This  result  deserves  attention. 

When  artificial  viscosity  is  not  considered,  the  stability  for 
a  one-dimensional  Lagrangian  coordinate  difference  grid  is  inde¬ 
pendent  of  the  cross-sectional  area  A(r),  and  the  Courant  condi¬ 
tion  still  holds,  i.e.,  A*<Ar/«[3)  .  But  since  the  artificial 

viscosity  term  is  added,  in  actual  computation,  we  multiply  it  by 
the  coefficient  0.6,  i.e.,  ai  — min  0.6  ^4. 

•i  ’ 

2.  NUMERICAL  RESULTS  AND  ANALYSIS 

Figure  2  shows  the  distribution  of  the  parameters  when  the 
major  shock  is  situated  in  the  middle  of  the  changing  cross-section. 
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Figure  2.  Axial  distribution  of  flow  parameters  (section 
with  changing  cross-section) 

Key  : /--medium  pressure  section;  2  —  contact  surface;  3  —  section 

with  changing  cross-section;  4 — new  shock;  5 — major  shock; 
6 — experimental  section;  7 — pressure;  8 — velocity;  9 — 


Figure  3.  Axial  distribution  of  the  flow  parameters 
(experimental  section) 


Key:  1 — medium  pressure  section;  2 — section  with  changing  cross- 
sectional  area;  3 — new  shock;  4 — experimental  section; 

5 — contact  surface;  6 — major  shock 


New  shock  wave  appears  behind  the  major  shock.  The  low  pressure 
and  low  density  regions  form  before  the  new  shock  (situated  to 
the  left  of  the  major  shock)  because  of  the  rapid  expansion  of 
the  gas  flow  in  the  conical  tube  when  the  separating  surface  has 
not  yet  entered  the  section  with  changing  cross-section  although 
the  major  shock  is  already  in  the  middle  of  that  section.  The  flow 
diagram  at  this  time  agrees  with  the  analysis  in  [1]. 
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Figure  4.  Computational  result  for  the  2  gas  models 
(pH=300kg/cm2 ) 

1 — medium  pressure  section;  2 — section  with  changing  cross- 
sectional  area;  3 — experimental  section;  4 — contact  surface 
5 — simplified  ideal  gas  model;  6 — variable  specific  heat 
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Figure  5.  Comparison  of  calculated  peak  shock  pressure 
with  measured  results 


1 — experimental  measured  value;  2 — measured  value  A;  3 — 
measured  value  0;  4 — measured  value  0;  5 — measured  value  + 
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Figure  3  is  the  situation  when  the  major  shock,  new  shock 
and  the  separating  surface  of  the  pauses  all  are  close  to  the 
experimental  section.  The  motion  of  the  separating  surface  of  the 
two  gases  has  already  overtaken  the  new  shock  when  t  =  7.75. 

From  the  diagram  one  can  see  that  to  accurately  calculate 
the  complete  flow  field  at  the  experimental  position,  it  is  necess¬ 
ary  to  consider  the  two-gas  media  model  in  front  of  the  changing 
cross-section.  This  has  not  been  discussed  in  [1].  In  the  work 
of  Chisnell  et  aL,  the  new  shock  and  the  separating  surface  between 
the  two  gases  are  ruled  out.  Therefore,  it  is  also  not  applicalbe 
for  this  computation. 

Figure  4  are  the  two  pressure  curves  calculated  with  the  two 

state  equation  schemes  under  the  condition  that  the  initial  explo- 

2 

sive  gas  excess  pressure  is  300  kg/cm  .  The  two  curves  basically 
overlap. 

Figure  5  is  the  comparison  between  the  calculated  result  and 
experiment.  The  calculated  result  of  the  peak  shock  pressure 
generally  differs  from  the  experimental  result  within  10%. 

Figure  6  is  the  curve  of  the  shock  tube  pressure  that  decreases 

with  time.  The  pressure  is  calculated  at  four  fixed  positions 

r  —  37.9.  55.0,  57.8,  62.5  when  the  initial  explosive  gas  pressure  is 

2 

100  kg/cm  .  It  basically  agrees  (Figure  7)  with  the  measured 
pressure  curve. 

Conservation  of  total  mass  and  total  energy  has  been  used  to 
test  the  accuracy  of  the  calculation  in  the  whole  computational  pro¬ 
cess.  When  the  dimensionless  time  t  =  4.66,  the  corresponding 

change  in  the  total  mass  is  0.34x10  ,  and  the  corresponding  energy 

_2 

change  is  -0.9x10 
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